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1.  INTRODUCTION 


This  report  presents  the  theoretical  background  and  computer 
programs  for  the  analysis  techniques  that  were  used  in  evaluating  data 
obtained  in  the  electromagnetic  pulse  (EMP)  testing  of  military 
communication  and  weapons  systems  under  the  PREMPT  program. ^ 

The  data  presented  are  initially  obtained  as  a voltage-versus-time 
trace  photographed  on  Polaroid  film.  This  trace  is  then  digitized  and  a 
time  series  of  digital  values  is  produced.  The  data  are  then  processed 
in  a digital  computer.  The  various  techniques  employed  in  reducing  and 
transforming  the  data  are  grouped  under  the  generic  title  "signal 
analysis."  Section  2 of  this  report  gives  a detailed  description  of  all 
the  algorithms;  also,  it  contains  complete  instructions  on  how  to  use 
the  signal-ainalysis  program. 

Another  technique  used  in  the  data-reduction  process  is  to 
represent  the  EMP  waveform  by  a set  of  parametrized  functions.  This 
technique  involves  a least-squares  fitting  procedure  which  is  discussed 
in  section  3.  Also,  section  3 contains  complete  instructions  on  the  use 
of  a least-squares  fitting  program. 

This  report  is  not  intended  to  be  exhaustive  on  the  subject  of 
signal  analysis  but  rather  to  present  to  the  EMP  community  a basic 
software  package  that  will:  (1)  accomplish  most  of  the  data  reduction 
for  EMP  work  and  (2)  be  easily  modified  to  include  any  additional 
techniques. 

2.  SIGNAL  ANALYSIS 

This  section  presents  the  theoretical  background  and  computer 
implementation  of  a number  of  techniques  for  reducing  and  transforming 
digital  time  series  produced  in  EMP  tests  under  the  PREMPT  program.  A 
complete  program  listing  annotated  with  comments  is  given  in  this 
section.  Several  versions  of  this  program  have  been  implemented  on  both 
the  IBM  370-195  and  CDC  6500  computers.  All  of  the  programming  was  done 
in  FORTRAN.  Some  of  the  subroutines  have  been  coded  in  assembly 
language  for  the  NOVA  minicomputer  but  are  not  reported  here.^ 


I ^rhe  PREMPT  program  Is  a joint  NMCSSC/DNA  effort  to  determine  the 
response  of  DCS  to  electromagnetic  pulses  generated  by  a high-altitude 
^nuclear  burst. 

I ^Further  details  are  presented  in  ''The  Interactive  Digitization  and 
^Editing  System  (IDES)''  by  Dr.  Thomas  A.  Tumolillo,  USA  Harry  Diamond 
Laboratories/  Adelphi/  MD  20783,  (Aug  1973) . 
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2.1  Data  Preprocessing 


Under  the  generic  title  "Data  Preprocessing"  is  included  the  many 
minute  details  that  are  necessary  to  prepare  the  raw  input  data  as 
obtained  from  a digitization  of  the  waveform,  so  that  it  is  suitable  for 
transformation  to  the  frequency  domain. 

The  following  sections  discuss  the  method  of  referencing  the  trace 
to  the  scope  graticule,  scaling  the  data,  time  ordering,  bit  reversal, 
and  a few  of  the  simpler  interpolation  schemes.  The  software  for  two  of 
the  simpler  interpolators,  the  linear  Lagrange,  and  the  linear  least 
squares  are  presented.  Higher  order  interpolators  have  been  used 
occasionally  in  the  PREMPT  program,  but  are  not  included  here  because 
the  simpler  methods  usually  work.  Simililarly,  no  techniques  for  time 
tying  of  digital  records  are  given.  Only  the  software  are  presented  for 
the  most  commonly  used  grid  and  tablet  referencing  schemes,  even  though 
some  of  the  more  complicated  procedures  are  discussed. 

2.1.1  Referencing  the  Trace  to  Graticule  and  Digitization  Tablet 

One  important  function  of  the  data  reduction  is  the  determination 
of  the  rotation  angle  of  the  scope  graticule  with  respect  to  the 
digitization  tablet,  the  zero  point  (origin)  of  the  graticule,  and  the 
scale  factors  for  the  X and  Y axes  in  the  graticule  (or  grid)  coordinate 
system.  There  are  many  possible  schemes  that  can  be  used  to  determine 
these  factors;  the  most  general  method  will  be  discussed  in  this 
section.  The  software  exists  for  the  general  procedure  as  well  as  for 
the  simpler  specific  case  implemented  in  this  signal  analysis  package. 

For  simplicity,  but  with  no  loss  of  generality,  let  the  grid 
points  be  symmetrical  about  the  origin  of  the  grid  coordinate  system 
(fig.  1);  call  this  origin  (X  ,Y  ).  In  the  grid  coordinate  system  the 
measure  grid  points  are  designated  by  the  arrays  XG(I)  and  YG(I).  In 
the  tablet  coordinate  system  the  measured  grid  points  are  designated  by 
the  arrays  XT  (I)  and  YT(I).  The  coordinate  systems  are  shown  in 
figure  1. 

The  transformation  equation  between  the  two  systems  is 


COS0 

sin6 


-sin0| 

cosoj 


(1) 


B 


Yt  Yg 


Figure  1.  Diagram  of  coordinate  systems. 

Applying  this  transformation  to  the  measured  grid  points,  we  have 
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The  points  are  symmetrical  in  the  grid  system;  thus,  we  have 
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N 

1 — 

X 

o 

II 

E 

1=1 

N 

XT  (I) 

Y 

o 

E 

1=1 

YT  (I) 

_ _ 



(3) 


It  is  convenient  to  shift  the  origin  of  the  tablet  coordinate  system  to 
the  ' thus, 


For  a specific  set  of  points 
transformation  can  be  linearized, 
we  have 


jcos0 

-sin0 

1 

^G 

|sin0 

COS0| 

1 

in  the  grid  system,  the  above 

Take  the  line  defined  by  Y = 0;  then 

G 


= cose  = 0 


(5) 


Let  a = COS0,.  then  X^  = best  value  of  0 can  be  estimated  from 
the  measured  points  along  the  grid  X axis  by  using  a least-squares 
technique-  Let 


- a XG(I)  ; 
I 


(6) 


minimizing  with  respect  to  a,  we  have 


a(x') 

9 a 


^XT(I)*XG(I) 

0 =i>  a = COS0  = 

^XG(I)  *XG(I) 
I 


(7) 


The ^ means  that  we  only  sum  over  those  points  on  the  grid  X axis. 

We  could  now  suitably  define  other  straight  lines  on  subsets  of  the 
measured  grid  points  and  get  further  estimates  of  0.  It  is  more 
convenient  to  use  a nonlinear  least-squares  technique,  and  extract  the 
best  value  of  8,  by  using  all  the  points  at  once. 
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(8) 


Let  Z(I)  = 


XT  (I) 
YT(I). 


•{■(O,!)  = 


(COS0  -sin0) 


(sin0  COS0 


and 


poun  ^ 

Lyg(i)J 


(9) 


X2  = ^ |Z(I)  - 1(0,1) |2. 

1=1 


(10) 


We  minimize  with  respect  to  0 by  requiring  9(x^)/98 
we  have 


= 0;  thus, 


= t 


jJxll  = 

30 


1=1  ‘* 


3?''' (9, 1) 
30 


(Z(I)  -<t>(0,D)  + (z(i) 
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1* 

Here,  is  the  adjoint  of  <*>. 


Assume  now  that  have  defined  an  iterative  process  for  evaluating 

6 and'l>,  at  the  iteration,  we  assume  that  is  given  by 


$(9,I)'v^  = ? + 


,_k 


39 


. A0 


48^  . - 6*= 


(12) 


By  substituting  equation  (12)  into  equation  (11),  a recursion  relation 
is  obtained  for  0.  It  can  be  shown  after  a slight  algebraic 
manipulation  that 


G^*  - C^sinO*^  + C 


:os0 


(13) 
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where 


n 

X!  (XG(I)*XT(I)  + YG(I)*YT(D) 


(XG(I)*XG(I)  + YG(1)*YG(D) 

1=1 


(14) 


N 

^ (XG(I)*YT(I)  - YG(I)*XT(I) ) 

S “ 1^1 (15) 

N ' 

Y {XG(I)*XG(D  + YG(I)*YG(D) 

1=1 


The  recursion  relation  generally  converges  quickly^  as  long  as  the 
initial  estimate  for  0 is  close  to  the  true  value.  We  used  fc^r  0^  the 
value  determined  in  equation  (7).  The  program  considers  tiiat  the 

iterative  scheme  has  converged  as  long  as  A0^  < 10  ^ and 


The  scale  factors  for  the  X and  Y axes  can  be  easily  calculated. 
For  illustrative  purposes,  suppose  we  measured  28  grid  points  as  shown 
in  figure  2. 


Initially,  the  arrays  XG  and  YG  are  defined  in  a DATA  statement  so  that 
the  points  IXG(I),  YG(I)1  corresx>ond  to  a grid  so  that 

XG(I)  r.  {-1.  . . . ,+.8,+l.  } 

YG(I)  L {-l.,-.75,...,4-.75,-fl.} 

so  that  for  example  [XG(27),  YG{27)]  = (+.4, -.5).  In  the  program,  one 
of  the  initial  redefinitions  of  the  arrays  we  make  is 

XG(I)  > XS*XG(I) 

YG(I)  ^ YS*YG(I) 
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Scale  factors  XS  and  YS  must  be  determined  (fig.  3), 
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Figure  3.  Rotation  between  G and  T systems. 

For  the  [joints  on  die  axis  (I  = we  have 

u 

[XS*XG(1)]^  = XT(I)*XT(I)  + YT(I)*YT(I)  . 
The  avera<je  value  for  XS  is 
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Figure  2.  Measured  grid  points. 
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(16) 


Similarly,  from  the  points  on  the  Y^  axis  (I  = 12,  20) , we  have 


f 


20 

Y,  (XT(I)*XT(I)  + YT(I)*YT(D) 

YS  = 1^12 (17) 

20 

Y YG(I)*YG(I) 

\|  1=12 


A complete  grid  measurement  is  rarely  carried  out  as  a routine 
operation  in  the  reduction  of  EMP  data.  Occasionally,  it  will  be  done 
to  test  the  linearity  of  the  system.  The  most  common  method  is  to 
measure  two  points  on  each  axis  and  calculate  all  quantities  from  these 
numbers.  The  following  equations  wore  implemented  in  the  signal 
analysis  package.  The  routine  calculates  the  grid  rotation  angle,  the 
scale  factor,  and  the  origin  point  from  two  measured  points  on  both  the 
X and  Y grid  axes.  Initially,  the  program  reads  in  the  coordinates  of 
the  two  points  measured  on  the  X axis — XI  and  X2 — and  on  the  Y axis — Yl 
and  Y2.  Then  the  four  grid  points  are  read  in  the  same  order  that  the 
coordinate  locations  were  read  in.  The  four  grid  points  are  stored  in 
the  arrays  XT (I),  YT(I),  I = 1,4.  The  origin  of  the  grid  (XO,  YO)  is 
given  by 


XO  = |X2*XT(1)-X1*XT(2) I /|xl-X2|  , 

YO  = |Y2*YT(3)-Y1*YT(4) j /1y1-Y2|  , (18) 


The  scale  factors  XS  and  YS  are  given  by 


XS  =/(XT(2)-XT(l) )**2  + (YT(2)-YT(1) ) **2  /(X2-X1)  , 

YS  =/(XT(4)-XT(3) )**2  + (YT (4 ) -YT (3 ) ) * *2  / (Y2-Y1) , 


(19) 
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The  rotation  angle  is  determined  by 


. _ n c*  /yT(2)-YT(l)  XT(4)-XT(3)\ 

tan(0)  0.5  ^XT(2)-XT(1)  " YT (4) -YT (3) / ' 

CT  = cos(e)  = l.//l+tan(0)**2  , 

ST  = sin(0)  = cos (0) *tan (0)  . (20) 


Then  the  subroutine  reads  scale  factors  T (nanosecond/grid  division)  and 
V (volts/grid  division),  then  recomputes  -the  scale  factors  as 


T^/XS,  VV/XS  , (21) 

which  have  units  nanoseconds/tablet  counts  and  volts/tablet  counts. 
Then  the  (x,y)  coordinates  of  the  trace  are  transformed  as. 


X(I)  - (CT*  (X(I) -XO)  4-  ST*  (Y(I)-Y0)-X2;)*T  , 

Y(I)  - (-ST* (X(I)-XO)  + CT* (Y(I)-Y0)-YZ)*V  , (22) 


where  (XZ,  YZ)  are  the  rotated  coordinates  of  the  "zero  point  of  the 
trace,"  that  is,  the  point  on  the  trace  at  which  the  signal  starts, 


XZ  ^ CT*(XZ-X0)  + ST*(YZ-Y0) 

YZ  - -ST*(XZ-X0)  + CT* (YZ-YO)  (23) 


2.1.2  Time  Ordering,  Bit  Reversal,  and  Interpolation  Schemes 

Time  ordering  of  the  array  is  a necessary  procedure  in  order  to 
remove  errors  introduced  in  the  digitization  x^rocess.  Occasionally, 
there  will  be  errors  in  the  grid  measurements  that  cause  portions  of  the 
trace  to  fold  back  in  the  time  sense  after  it  is  rotated.  Similarly, 
inaccurate  movement  of  the  digitization  operator's  hand  may  also  cause  a 
few  points  to  be  folded  back  in  the  time  sense.  In  some  cases  the 
digitization  hardware  will  allow  consecutive  digital  points  to  have  the 
same  time  value.  If  these  measurement  ambiguities  are  not  removed,  they 
will  cause  considerable  error  in  the  high-frequency  part  of  the 
transforms.  This  correction  of  the  data  is  handled  in  subroutine  CST 
OUT.  This  routine  casts  out  those  points  in  the  array  that  are  folded 


15 


back — that  is,  if  XF  is  the  input  time  array  and  XF(K)<XF(I)  for  K = I + 
1,  1 + 2,...  then  SF(K)  is  deleted  from  the  array.  Similarly,  if  XF  (K) 
= XF(I)  for  some  set  of  K then  the  routine  averages  the  amplitude  YF (K) 
to  create  a single  value  at  that  value  XF(I). 

Bit  reversal  of  the  array  refers  to  a specific  reordering  of  the 
elements  of  a digital  time  series  prior  to  its  entering  the  FFT 
subroutine.  It  is  done  so  that  after  transformation  the  frequency 
domain  arrays  are  in  ascending  order  of  the  frequency  value.  The  term 
bit  reversal  arises  from  representing  the  index  of  an  array  I in 
base-two  notation.  For  example,  suppose  we  have  the  65th  element  of  a 
1024  element  array,  then  65^^  = OOOIOOOOOI^.  The  reverse  of  the  number 

is  1000001000^  = 520^q.  To  bit  reverse,  we  swap  the  elements  65  and  520 

of  the  original  array. 

In  subroutine  LNYQ,  the  Lagrangian  methods  for  interpolation  are 
used.  A brief  description  of  Lagrange  interpolation  is  given  here. 

It  is  generally  assumed  that  the  function,  f,  interpolated  here 
behaves  like  a polynomial;  thus,  in  order  to  calculate  f approximately 
at  a point  x,  we  find  a polynomial  approximation  g for  f good  in  the 
neighborhood  of  x.  Lagrange  showed  that  there  is  a unique  polynomial  of 
degree  n having  n + 1 values  f^  at  n + 1 distinct  points  x^  f i = 

0,.«.n.  That  polynomial  is  g^, 

n 

n • (24) 

i=o 

j=o 

j^i 


For  the  software  presented  in  this  package,  n is  restricted  to  the  value 
1.  Thus, 


g^  (x) 


f(x  )-f(fx  ) 

1 o_ 


X + 


H Cl-x  + C2  . 


(25) 


In  the  program  the  function  f is  called  YF(I),  x.  is  replaced  by  the 
time  array  XF(I),  and  the  interpolated  values  are  put  into  the  real  part 
of  the  complex  array  YNYQA (K) 
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ri-:al  [ynyoa(k)1  = ci-x  + c2, 


YF  (I)-YF(I-l) 
XF(I)-XF(I-1)  ' 


C2 


XF(I) *YF(I-1)"XF  (I-I) *YF(1) 
XF(I)-XF(I-1) 


(26^ 


In  subroutine  NYQST,  a linear  x^^>lynomial  is  fitted  to  a set  of 
points  in  the  arrays  YF(I),  XF(I),  I = LB,...,LT,  After  determining  the 
polynomial,  the  program  evaluates  it  at  the  ^^redetermined  interpolation 
point. 


The  theory  behind  the  least  squares,  polynomial  fitting  x)rograms 
is  straightforward.  We  need  to  minimize  with  respect  to  the  C(J) 
where 


X 


2 


LT 

Y,  Iyf(i) 


I=LB 


M 

~Y  C(J)*XF(I)**J]2  , 

J=0 


(27) 


Setting 


dim 

0C(J) 


0 J = 0, ,M  yields 


B = A-C 


(28) 


where  the 


(K,L) 


th 


element  of  the  matrix  A 


is 


LT 

A(K,L)  = Y XF(I)**(K  + L - 2),  K,  L = (29) 

I=LB 


th 

the  K element  of  the  vector  B is 
LT 

B(K)  = Y YF(I)*(XF(I)**(K  - D)  K=1,...,M,  (30) 

I=LB 

and  the  element  of  the  vector  C is  just  the  polynomial 

coefficient.  Upon  inversion  we  find  the  solution  for  C 


C = 


B 


(31) 
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Both  of  the  interpolator  subroutines  take  due  account  of  the  end 
points  of  the  arrays  and  minimize  the  number  of  calculations  when  more 
than  one  interpolated  point  falls  between  the  same  time  values. 

2.2  Digital  Filtering 

This  section  briefly  reviews  the  theory  of  digital  filters. 
Frequently  in  EMP  work,  the  signal  is  contaminated  by  high-frequency 
noise  arising  both  from  the  nature  of  the  measurements  and  the 
digitization  process  which  reduces  the  continuous  signal  to  a digital 
record.  This  high-frequency  component  can  generally  be  eliminated  by 
passing  the  digital  record  through  a low-pass  digital  filter.  Another 
useful  application  of  low-pass  filters  is  in  reducing  the  number  of 
digital  values  needed  to  accurately  calculate  the  Fourier  transform  of  a 
waveform  at  low  frequencies.  For  example,  if  there  is  a signal  with 
frequency  content  up  to  250  MHz,  the  Nyquist  criterion  is  satisfied  by 
sampling  the  signal  every  2.0  nsec.  If  the  signal  has  a 2-usec 
duration,  then  1000  numbers  must  be  stored  for  the  Fourier  transform 
routine.  However,  if  most  of  the  significant  frequency  content  is 
contained  in  a frequency  band  up  to  say  50  MHz,  then  after  filtering, 
only  200  numbers  must  be  stored  to  adequately  represent  the  signal  and 
obtain  the  Fourier  transform  without  worrying  about  foldover  effects. 
Bandpass  and  high-pass  digital  filters  have  application  in  EMP  work  when 
one  is  interested  in  studying  only  a certain  region  of  the  frequency 
spectrum  so  that  correlations  between  equipment  upset  and  damage  and  the 
induced  signal  can  be  determined. 

2.2.1  Theory  of  Digital  Filtering  and  Z Transforms 

In  EMP  work,  there  is  generally  (after  interpolation  or  as  a 
result  of  the  digitization  process)  a sequence  of  numbers  u(k),  k=0,...N 
that  must  pass  through  a linear  discrete  system  (generally  a difference 
equation)  in  order  to  limit,  in  some  manner,  the  frequency  content  of 
the  signal.  The  output  of  the  linear  system  is  denoted  here  by  y(k), 
k=0,...,No  By  a linear  discrete  system  is  meant  a system  in  which  the 
output  y(k)  is  expressed  as  a linear  combination  of  inputs  and  past 
outputs;  thus. 


n 

y(k)  t 

j=;i 


m 

b(£)u(k-il) 


(32) 


and  y(k),  u(k)  =0,  k < 0. 


IR 


The  Z transform  is  used  to  simplify  the  analysis  and  synthesis  of 
the  digital  filter  represented  by  equation  (32) — for  example,  generate 
the  set  of  constants  a(j)  and  b(X.). 

The  Z transform  of  a sequence  of  numbers  f(k),  k=0, — N f (k)  = 0, 
k < 0 is  defined  by 


Z[f(k)]  = F(z) 


Z f(k)z 
k=0 


(33) 


where  z is  an  arbitrary  complex  number. 

The  Z transform  of  the  input-output  signals  are  related  to  one 
another  by 


Y(z)  = H(z)  U(z) 


(34) 


where  U(z)  and  Y (z)  denote  the  Z transforms  of  the  input  and  output 
signals  respectively,  and  H(z)  is  the  system  transfer  function,  and  is 
given  by 


m 


H(z)  = t>(2.)z'^/ll  + a(j)z 

2,=0  \ j = l 


-1 


(35) 


-k 

Proof:  Multiply  each  side  of  equation  (32)  by  z and  sum  over  all  k. 


y(k)z  ^ a(j)  ^y(k-j)z 

k=0  3=1  k=0 


b(^)  ^u(k-f.)z”  . (36) 

Z=0  k=0 


Using  the  properties  of  the  Z transform,  this  can  be  rewritten  as 


n _ . m _ j 

Y (z)  + ^ a(j)z  ^ Y(z)  = ^ bC^Oz 
j = l 
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U(Z)  ; 


(37) 


or , as 


n 

1 + ^ a ( j ) z ^ 

' o' 

Y(z) 

= 

y b(i)z 

j=l 

o' 

\ 

(38) 


Thus,  equation  (35)  follows  and  may  be  rewritten  as 


H (z) 


b n (z-55. ) / n 

i=l  1=-! 


(39) 


Here,  z.  and  p.  are  called  the  zeros  and  poles  of  the  system  transfer 
function  H(z). 

One  of  the  most  important  properties  of  the  transfer  function  is 
the  fact  that  the  location  of  the  zeros  and  poles  has  an  enormous  effect 
on  how  the  system  transmits  different  types  of  inputs.  Thus,  a system 
can  be  synthesized  that  will  pass  some  inputs  and  reject  others  by  a 
X)roper  selection  of  the  zeros  and  poles. 

Assume  the  application  of  a sinusoidal  input  u(k)  = sin(kwT)  k=0, 
1,  2,...  to  our  systems.  Then  the  resultant  steady-state  response, 
y(k),  is  given  by 


y(k) 


sin  (kv:T-»-0 ) 


(40) 


Here,  the  sinusoidal  response  of  the  system  is  obtained  by  evaluating 

iwT 

the  system  transfer  function  H(z)  at  z=e  , where  w is  the  radian 
frequency  of  the  input  sinusoid  and  T is  the  underlying  sampling  period. 

0 is  the  phase  angle  of  H (e 

prove  equation  (40)  let  u(k)  = sin(kwT)  k = 0,  1,  2, 

given  by 


\ . 1 iwT\ 

/ iwT  \ 

i ^ 

) — that  IS,  H(e  / 

- 

H(e  ) 

o 

To 

then,  U(z)  is 


U (z)  = 


zsinwT 


(41) 
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The  Z transform  of  the  system  response  is  given  by, 


Y(z) 


H (z) zsinwT 


(4  2) 


Since  only  stable  systems  are  considered,  all  the  poles  of  H(z) 
must  be  inside  the  unit  circle.  Therefore,  none  of  the  poles  of  H(z)  is 
iv^T  ~iwT 

at  e or  e .A  partial  fraction  expansion  of  equation  (42)  yields 


Y(z)  = 


az  + bz + terms  due  to  poles  of  H(z),  (43) 

iwT 


\z-e 


where 


a = and  b = ~H(e  ^'^^)/2i  . 

Noting  that  = n(e  and  writing  = Me^*^',  we  find  that 


M 

Y(z)  = - 


ze 


i0 


~ze 


-i0 


terms  due  to  poles  of  H(z)  . (44) 


Taking  the  inverse  Z transform  of  equation  (44)  we  obtain 


y(k) 


M f i(kwT-H0) 
2i 


-i (kwT+0) 

e J + transient  response 

generated  by  poles 
of  H(z)  ; 


(45) 


or. 


y (k)  = M sin(kwT+0) 


+ . ^(k) 

transient 


(46) 
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In  the  steady-state  y transient  (k)  - 0 as  k becomes  large.  Thus 


y(k) 


sin (kwT+0) 


(47) 


The  quantity  M=|H(e^'^  )|  is  called  the  system  gain  factor.  To 
filter  out  a given  sinusoid,  pick  a system  -so  that  M'vQ;  or,  to  amplify  a 
given  sinusoid,  pick  M>1. 


As  has  been  observed  in  equation  (39) , the  transfer  function  is  a 
ratio  of  polynomials  in  the  variable  z.  Therefore,  the  system  gain 
factor  that  is  equal  to  the  magnitude  of  the  transfer  function  evaluated 
at  z=exp(iwT)  may  always  be  expressed  as  a ratio  of  polynomials  in  the 
variables  cos(wT)  and  sin(wT).  Thus,  different  filters  can  be 
synthesized  by  investigating  ratios  of  trigonometric  functions.  For 
example,  a low-pass  filter  with  half-power  point  w,  has  the  following 
squared  gain  factor. 


1 

^ ^ tan^"(wT/2) 
tan^'^|wjT/2) 


(48) 


By  a considerable  amount  of  algebraic  manipulation  equation  (48)  can  be 
written  as 


tan^^|w^T/2j  |l+zj^^ 


(49) 


wliere  z = exp(iwT),  and  the  2n  poles  p^  are  given  by 


1 - tan^l 

[w,T/2| 

1 -f  y-1  2tan|w^T/2\sinO^ 

1 - 2tan| 

|WiT/2^ 

1 cos0^  + tan2 1 

|w^T/2] 

1 ^ 

(50) 


where, 


0^  = (i-1)  71/2,  n odd 

= (2i-l)  tt/2,  n even  . 

It  can  be  shown  that  of  the  2n  poles,  p.,  exactly  n lie  inside  the 
unit  circle  and  n outside.  Let  p^ , ^ poles  inside 

the  unit  circle.  The  transfer  function  that  has  the  desired 
squared-gain  factor  is  given  by 


H(z) 


b (1-Hz)^ , 


<51) 


where  b is  chosen  so  that  the  steady-state,  unit-step  response  has 
magnitude  one,  H(l)=l;  thus. 


b 


(52) 


The  remaining  iK>les  ^n*f 2' ’ ‘ ' '^2n  with  |H(z)|^can  be 
shown  to  arise  from  the  process  of  determining  the  squared-gain  factor. 

From  the  foregoing,  a well-defined  procedure  exists  for 
synthesizing  a low-pass  filter.  It  can  be  summarized  in  the  following 
steps: 


(a) 

Determine  the  half-pcwer  point 

(b) 

Determine  the  value  n — using  equation 
gain  at  frequency  w^, 

(17)— by 

specifying  the 

(c) 

Find  the  n roots  p given  by 

satisfy  il^|<  1,  and 

equatior 

. (A9)  which 

(d) 

Determine  the  difference  equation  which  has 
function  given  by  equations  (51)  and  (52) . 

the  transfer 

A squared-gain  factor  that  corresponds  to  that  of  a high-pass 
filter  is  given  by 
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1 


(53) 


^ ^ cot'^'^(wT/2) 

cot^"  (W2T/2) 


Here,  the  half -power  x)oint  is  denoted  by  . Formulas  anaiagous  to 
equations  (50)  and  (51)  may  be  derived  Then,  it  the  poles  and  zeros  of 
a low-pass  filter  with  half-power  point  (tt/T-w^)  are  rotated  through 
TT  radians  in  the  complex  plane,  the  pole-zero  pattern  of  a high-pass 
filter  is  obtained  with  half-power  point  w^  and  the  same  gain-factor 
falloff  outside  its  passband  is  obtained.  Thus,  if  we  have  H(z)  for  the 
low-pass  filter  given  by 


n 


then  the  high-pass  filter  is  given  by  H*  (z) 


H’  iz) 


(54) 


(55) 


2.2.2  Examx:)les  of  Filtering  FMP  Data 

To  illustrate  the  implementation  of  the  algorithms  described  in 
section  1,  a typical  EMP  waveform  was  selected  from  the  vast  amount  of 
data  collected  at  the  Polk  City  AUTOVON  EMP  tests  and  processed.  Figure 
4 plots  the  digitized  waveform  after  it  has  been  digitized,  time 
ordered,  and  interpolated  at  1.63-nsec  intervals  using  a Lagrange 
interpolator. 

Figure  5 plots  the  power  spectrum  after  passing  the  digital  record 
through  a fast  Fourier  transform  routine.  All  the  power  is  contained  in 
two  peaks  at  13.2  and  24.0  MHz.  Several  filters  were  then  synthesized 
and  the  digital  record  passed  through  filters  before  processing  it 
through  the  fast  Fourier  transform  routines.  Table  I lists  the 
half-power  points,  w , the  gain  at  the  higher  frequency  w^  used  to 
determine  n,  and  tne  filter  coefficients  a(j),  j=l,...,n,  b(j), 
j=0.L.,n,  which  were  calculated  for  each  synthesized  filter.  Figures  6 
through  12  are  plots  of  the  power  spectra  obtained  by  using  a fast 
Fourier  transform  routine  on  the  f iltered-time  series.  The  general 
result  is  fairly  evident  from  these  plots — namely,  that  as  the  number  of 
poles  is  increased  and  consequently  the  gain  rolloff  at  the  half-power 
point  is  increased,  a sharper  filtering  is  obtained. 

Although  it  is  not  apparent  from  the  plots  of  the  power  spectra 
presented,  there  is  much  greater  definition  of  the  peaks  in  the  power 
spectrum  relative  to  the  noise  background.  The  high-frequency  noise  was 
reduced  by  a factor  of  100. 
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POWER 


FREQUENCY  (UNIT  1. 198  MHi) 


Figure  5.  Plot  of  power  spectrum  obtained  without  filtering. 
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TABLE  I.  PARAMETEPS  CALCULATED  FOR  EACH  SYNTHESIZED  FILTER 


Filter 

Ho. 

Half*power 
point  W 
(MHz)^ 

Gain  H- 

“2 

Sai^pling 
Interval  T n 
(p«ec) 

1 

20 

0.4 

25 

1.63  1 

2 

20 

0.2 

25 

1.63  4 

3 

20 

0.1 

25 

1.63  5 

4 20  0.07*  25  1.63  6 


5 20  0.05  25  *‘6.52  6 


X8  0.2  25  1.63  • 3 


7 17.5  0.05  25  1.63  5 


a(j) ,j-l,...,n  b(J) , j-1,. ..n+l 


-0.8135657E00 


-0.3465070E  01 
0.4533676E  01 
-0.2652045E  01 
0. 584804 7E  00 

-0.4337327E  01 
0.7563597E  01 
-0.6625226E  01 
'0.2913780E  01 
-0.5145599E  00 


-0.5208673E  01 
0.1135041E  02 
-0.1324103E  02 
0.8718896E  01 
-0.3071778E  01 
0.4523243E  00 


-0.2848268E  01 
0. 384500 3E  01 
-0.2959958E  01 
0.1351626E  01 
-0.3426920E  00 
0.3745000E-01 


-0.2631725E  01 
0.2328194E  01 
0.6912119E  00 


-0.4420094E  01 
0.7844978E  01 
-0.6986225E  01 
0.3120767E  01 
-0.5592871E  00 


+0.9321 71 3E00 
+0.9321713E00 

0.8540331E-04 

0.3416131E-03 

0.5124197E-03 

0.8540331E-04 

0.8235322E-05 

0.4117661E-04 

0.8235322E-04 

0.8235322E*04 

0.4117661E-04 

0.8235322E-05 

0.79 374 38E-06 
0.4762463E-05 
0.1190616E-04 
0.1587487E-04 
0.1190616E-04 
0.4762463E-05 
0.7937438E-06 

0.1299385E-02 

0.7796306E-02 

0.1949077E-01 

0.2598769E-01 

0.1949077E-01 

O.7796306E-02 

0.1299385E-02 

0.6571114E-03 
0. 19713 34E-02 
0.1971334E-02 
0.6571114E-03 

0.4385640E-05 

0.2192819E-04 

0.4385639E-04 

0.4385639E-04 

0.2192819B-04 

0.438S640E-05 
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•Figure  6. 


FREOUENCY  (UNIT  1.198  MH2) 


Plot  of  power  spectrum  obtained 
after  passing  the  digital 
•record  through  filter  No.  1 
(table  I) . 


Figure  7.  Plot  of  power  spectrum 

obtained  after  passing 
the  digital  record 
through  filter  No.  2 
(table  I) . 
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POWER 


Figure  8. 


Plot  of  power  spectrum  obtained 
after  passing  the  digital 
record  through  filter  No.  3 
(table  I) . 


FREQUENCY  (UNIT  1190  MHi) 


Figure  9.  Plot  of  power  spectrum 

obtained  after  passing 
the  digital  record 
through  filter  No.  4 
( table  I ) . 
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FREQUENCY  (UNIT  I 190  MHi) 


POWER 


( 


Figure  10.  Plot  of  power  spectrum  obtaine 
after  passing  the  digital 
record  through  filter  No,  5 
(table  I) . 
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FREQUENCY  (UNIT  I 196  MHi) 


Figure  11.  Plot  of  power  spectrum 

obtained  after  passing 
the  digital  record 
through  filter  No.  6 
(table  I) . 
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Figvire  12.  Plot  of  power  spectrim  obtained  after  passing  the  digital 
record  through  filter  No.  7 (table  I ) . 


2,3  Transform  Techniques 

Three  different  popular  transform  techniques  are  discussed  below. 
2.3.1  Fast  Fourier  Transform 


The  fast  Fourier  transform  is  by  far  the  method  most  preferred  for 
generating  the  Fourier  transform  of  a digital-time  series.  A short 
discussion  of  the  method  is  given  in  this  section. 


F(u)  of  a function  f(t)  is  defined  as 


F(o)) 


J f(t)e"^‘^^  dt 


(56) 


If  f(t)  is  nonzero  only  over  a finite  time  interval  T,  then  it  is  a good 
approximation  to  write  the  Fourier  transform  as  a discrete  sum. 
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N-1 

F(u>)  = At  X]  f (kAt)e"'^^'^^,  NAt 
k=o. 


= T . 


The  inverse  transform  is 


f (t) 


or 

= / 


F(w)  do) 


('i7) 

H* 


(•38) 


If  F(w)  has  no  (or  a negligible)  frequency  content  above  u)  the 

• - 1 ^ ^ max 

inverse  transform  can  be  written  as 


f (t) 


Aw 

2tt 


^ F(rAw)e  , NAw  = 2w 

^ max 

r=o 


(S9' 


For  a band- limited  signal,  we  know  from  the  Nyquist  criterion  that  the 
sampling  interval  At  must  be  chosen  so  that  At  = 27t/2w^  . Similarly, 
for  a time-bounded  signal,  the  frequency-sampling  interval  must  be 
chosen  so  that  Aw  = 2w/T, 


AtAw  < 


T ^ 
N T 


il 

N 


(60) 


Thus, 


F(rAw)  = At  ^ f(kAt)e  ^ i;=o,  1,..^,N  - 1 . (61) 

k=o 

For  notational  convenience,  we  write  F^  = F(rAw),  fj^  = f(kAt),  W = 
e“*2;Ti/N,  and  drop  the  factor  At. 


N-1 

'^r  = E 

k=o 


rk 


r = o, . . . , N - 1 . 


(62) 


Divide  the  time  series  of  points,  f^,  into  two  functions  and  h^^. 


N 

~ 0 / 1/  2 0 0 m 2 


1 


(6J) 


^2k' 


\ ^ ^2k+l, 


The  discrete  Fourier  transforms  of  g,  and  h,  are  G and  fl  , 

k k r r 

respectively .. 


(N/2)-l 


"r  = E 


k=0 


-2’'irk/(N/2) 


r = 0, 1,  — 2 - 1 


(G4) 


H = 
r 


(N/2)-l 

E 

k=o 


■2iiirk/(N/2) 


Equations  (61)  and  (62)  mav  be  rewritten  as 


(N/2)-l  -2Tirk/(M/2)  ^ -2-ir/N  <N/2)-l  -2r,irk/ (N/2) 

= E ^2k®  ^ E ^2k+l'^ 

k=0  k=0 


Usinq  equations  (63)  and  (64) , we  have 


F = G + V;  H 0 < r < N/2 

r r r “ 


^r+N/2  ” '^r 


W H 


0 <.  r < N/2 


(66) 


We  may  therefore  compute  the  Fourier  transform  of  a function  sampled  N 
times  by  evaluating  two  Fourier  transforms  of  the  function  sampled  N/2 
times.  The  computations  of  and  can  be  reduced  to  the  evaluation 
of  sequences  of  N/4  samples.  If  N = 2 r n such  reductions  can  be  made 
by  applying  equations  (63)  and  (66)  first  for  N,  then  N/2^  and  finally 
for  a two>-noint  function^  The  Fourier  transform  of  a one-point  function 


JJ 


is  just  the  sample  itself.  Beside  the  obvious  savings  in  computer 
running  time  by  using  this  successive  reduction  scheme/  the  transform 
can  also  be  done  in  place — that  is,  in  each  stage  of  reduction  the 
intermediate  results  are  written  over  the  original  array.  This 
in-place  reduction  requires  a rearrangement  of  the  original  array  called 
bit  reversal  (subroutine  SRTFUR) . A complete  rearrangement  and 
overwriting  sequence  for  an  eight-point  sampled  function  is  illustrated 
in  figure  13. 


Figure  13.  Rearrangement  and  overwriting  sequence  for  an  eight-point 
sampled  function. 


Each  arrow  in  the  diagram  means  that  the  term  at  the  origin  of  the  arrow 
must  be  added.  A variable  next  to  the  arrow  acts  as  a multiplier  to  the 
additive  term.  Thus,  at  the  first  overwriting,  sequence  f^  is  replaced 
by  is  replaced  by  ’*^f ^ . 

Subroutine  FFT  will  compute  the  transform  of  any  array  with  2^  elements 
as  described  above.  The  only  restriction  is  that  imposed  by  the  finite 
memory  size  of  the  computer  being  used.  Several  examples  of  the  power 
spectrum  derived  from  the  real  and  imaginary  parts  of  the  transform 
generated  by  FFT  are  shown  in  figures  5 through  12. 
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2o3.2  Fast  Walsh  Transform 


The  previous  section  transforms  the  digital-time  series  to  the 
frequency  domain  by  using  a particular  set  of  orthogonal 
functions — namely  sines  and  cosines « Another  set  of  orthogonal 
functions,  which  are  used  extensively  in  communications  theory,  are  the 
Walsh  functions  that  are  used  primarily  to  represent  logic  signals^, 
Their  most  appealing  feature  is  that  the  digital  Walsh  transform 
algorithm  is  about  an  order  of  magnitude  faster  than  the  Fourier 
transform  algorithmo 

The  Walsh  functions  are  wal(k,6),  sal(k,6),  and  cal(k,0). 


wal(2k,9)  = cal(k,0) 
wal(2k-l,0)  = sal(k,0) 


They  are  defined  on  the  time  interval  T,  9 is  the  normalized  time  0=t/T, 
and  k is  called  the  sequency.  The  sequency  is  equal  to  the  average 
number  of  zero  crossings  of  the  function  per  unit  time.  The  functions 
sal  and  cal  are  similar  to  the  sine  and  cosine  functions.  The  sequency 
of  the  Walsh  functions  plays  a similar  role  as  the  frequency  for  the 
sinusoidal  functionSo  One  definition  of  the  Walsh  functions  is  through 
a difference  equation. 


wal(2k+p,0)  = ^ ^ wal 

[f.,2(e  i) 

- ( 

- j^^^^walTk 

.2(e-  i) 

L 

L 

where 


k — 0,  1,  2,vw», 

[k/2]  is  the  largest  integer  less  than  or  equal  to  k/2, 
p = 0 or  1,  and 
wal(0,6)  = 1 0 1 j 

0 0 > i (68) 

Z V 


A few  of  the  Walsh  functions  are  shown  in  figure  14. 


Figure  14. 


Examples  of  Walsh 
imply  +1,  light 


functions  for  interval 
areas  imply  -1) . 


(dark  areas 
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A signal  f(t)  may  be  expanded  in  a Walsh  series 
f(t)  = a(o)wal(o,t)  + ^ [a  (k)cal(k,t)  + a (k)sal(k,t)l  , 

k=l  L ^ ^ J 

,-T/2 

a(o)  = J f(t)dt  , 

-T/2 

rT/2 

ac(k)  = J f(t)cal(k,t)  dt  , 

-T/2 

T/2 

as(k)  = j f (t)sal (k,t)dt  . (6g) 

-T/2 

Just  as  in  a Fourier  expansion,  the  sum  of  the  squares  of  the 
expansion  coefficients  give  the  sequency  energy  spectrum.  The  Walsh 
power  is 

E(k)  = a^  (k)  + (k)  . (70) 

To  evaluate  the  coefficients,  a fast  Walsh  transform  algorithm  can 
be  derived  that  is  similar  to  the  fast  Fourier  technique.  The  main 
difference  is  that  the  reduction  cannot  be  done  in  place.  The  steps  for 
an  eight-point-sampled  function  is  shown  in  figure  15.  The  arrows  have 
the  same  meaning  as  in  the  overwriting  sequence  for  the  FFT.  Thus,  at 
the  first  overwrite,  f^  is  replaced  by  f^+f^. 

For  the  EMP  data  of  figure  4,  the  Walsh  power  was  calculated  and 
the  results  plotted  in  figure  16«  The  same  figure  plots  the  Walsh  power 
after  filtering  the  data  through  filter  No.  7,  table  I. 


iV 


X 


Y~-^X 


X 


Qj(4) 


Figure  15.  Overwriting  sequence  for  Walsh  transfonn. 


Figure  16.  Walsh  power  for  the  time  series  of  figure  4. 


2,3.3  Refined  Spectral  Densities  and  Autocorrelation  Function 


In  EMP  work,  there  is  a great  need  for  implementing  numerical 
algorithms  for  the  processing  and  storage  of  digitized  waveforms,  which 
are  the  main  output  of  all  EMP  tests.  This  subsection  describes  one 
algorithm  to  generate  the  power  spectrum  of  a digital  record  from  its 
autocorrelation  function.  The  theory  presented  below  is  illustrated  by 
several  exeimples. 

A procedure  will  now  be  defined  for  estimating  the  power  spectrum  of  a 
uniformly  spaced,  discrete  time  series  of  finite  length. 

If  C(T)  is  the  autocovariance  function  for  a time  waveform  X(t), 
then  by  definition 

C(T)  = lim  - I X(t)X(t+T)dt  . (71) 

-T/2 

The  power  spectrum  P(w)  of  the  time  waveform  is  then  given  by 

oo 

P(w)  = 2 f cos(wt)  C(t)dt  . (72) 

0 

For  a uniformly  spaced  discrete-time  series  of  finite  length, 
denoted  by  X^,  X^,...,  X compute  the  mean  lagged  products,  C , with  lag 
interval  Ax=  hAt,  andAt  is  the  time  interval  between  adjacent  values  of 
the  time  series. 


C 


r 


1 

n-hr 


E 


q=0 


q qthr 


r=0,  1,. 


, m , m 


(73) 


Next,  compute  the  '*raw  spectral  density  estimates”  Vr. 


V 


r 


Ax 


C + 2 
o 


E 


q=l 


C cos  qrTT  + C cos 
q — — m 

m 


rtr 


(74) 


The  frequency  corresponding  to  r is  r/2mAx. 


to 


We  next  calculate  the  refined  spectral  density  estimates  according 


U = 0.23  V , + 0.54  V + 0.23  V . (75) 
r r-1  r r+1 

To  illustrate  the  implementation  of  the  algorithms  given  by 
equations  (73),  (74),  and  (75),  refer  once  again  to  figures  4 and  5. 
Figure  4 plots  the  digitized  waveform  after  it  has  been  digitized,  time 
ordered,  and  interpolated  at  1.63-nsec  intervals  with  a Lagrange 
interpolator.  Figure  5 plots  the  power  spectrum  after  passing  the 
digital  record  through  a fast  Fourier  transform  routine.  Most  of  the 
power  is  contained  in  the  two  peaks  at  13.2  and  24.0  MHZo 

Figures  17  through  23  are  plots  of  the  mean  lagged  products  and 
spectral  densities  for  different  values  of  the  lag  interval.  At.  It  is 
seen  that  all  of  the  frequency  content  of  the  power  spectrum  is 
accurately  calculated  until  the  lag  interval  exceeds  the  Nyquist 
sampling  rate,  A > t = l/2f^^,  where  f^  is  the  largest  expected 
frequency  content  or  the  record.  For  tiiese  data,  At^^  = 20  nsec. 
Figure  22  shows  the  power  spectrum  clearly  broadened  and  thus  fold-over 
effects  on  the  lower  frequency  peak.  Since  the  autocorrelation  function 
is  close  to  zero  for  times  greater  than  0.5  jjsec,  the  power  spectrum  can 
be  accurately  calculated  with  a lag  of  1.0  and  a lesser  number  of  mean 
lagged  products  (fig.  23). 
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POWER 


Refined  spectral  density 
estimates  for  the 
autocorrelation  func- 
tion of  figure  17. 


FREQUENCY  (UNIT  0.599  MH2) 


Figure  19.  Plot  of  mean  lagged 
products  C , r=l. 
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POWE 


FPEOUENCY  'UNIT  0 599  MHz) 


Figure  21.  Plot  of  mean  lagged 
products,  r=l 
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FREQUENCY  (UNIT  0 599  Mh/) 


QC 

UJ 


Figure  23.  Refined  spectral  density 
estimated , obtained 
after  truncation  of 
the  autocorrelation 
function  of  figure  6 
at  0.1502  ysec. 


Figure  22.  Refined  spectral  density 
estimates  for  the 
autocorrelation  func- 
tion of  figure  21. 
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3.  PROGRAM  DESCRIPTION  AND  INPUT  FORMATS 


A flow  chart  of  the  main  program  is  presented  in  figure  24.  Each 
subroutine  performs  a definite  signal  analysis  function.  As  indicated, 
the  user  is  privileged  to  employ  the  technique  of  his  choice  by  setting 
the  control  parameters  on  the  input  cards.  The  subroutines  and  their 
functions  are  listed. 


3.1  Subroutine  Descriptions 


MAIN:  This  routine  directly  or  indirectly  calls  the  rest  of  the 

subroutines  and  thus  controls  the  passage  of  the  program 
through  all  the  signal  analysis  options. 

READIN:  This  routine  reads  in  all  input  parameters  and  data. 

SCARTP:  This  routine  calculates  the  rotation  angle,  the  sine  and 

cosine  of  the  rotation  angle,  the  origin  (center)  of  the 
scope  graticule,  and  the  scale  factors  (tablet  units/div) 
along  both  the  time  and  voltage  axes. 

ROT:  This  routine  rotates  and  scales  the  input  time  and  amplitude 

arrays  given  in  digitization-tablet  units  into  units  of  time 
and  volts. 

CSTOUT:  This  routine  checks  the  time  ordering  of  the  digital  time 

series.  It  discards  those  points  from  the  array  whose  time 
values,  t^,  i = 1,...,N,  do  not  satisfy  > t^.  On  those 

points  that  have  the  same  time  values,  the  program  averages 
the  amplitude  values. 

LNYQ:  This  routine  interpolates  at  the  Nyquist  sampling  intervals 

with  a linear  interpolator.  If  the  time  value,  X,  satisfies 
XF(I-l)  < X < XF(I),  where  XF  is  the  time  array,  then  the 
interpolated  amplitude  is  Y = Cl  * X C2  where 

Cl  = (YF(I)-YF(I-1))/(XF(I)-XF(I-D) 

C2  = (YF(I-1)*XF(I)-YF(I)*XF(I-1))/(XF(I)-XF(I-1)  ) 


NYQST:  This  routine  interpolates  at  the  Nyquist  sampling  intervals 
by  using  a least-squares  interpolator o If  the  time  value, 

X,  satisfies  XF(I-l)  < X < XF(I)  where  XF  is  the  time  array, 
then  a least-squares  fit  of  the  function  $ = ajX(I)+a^  is 
made  to  the  set  of  points  (XF(J),  YF(J)),  J=I-2,  I-l,  I, 
1+1^  Thus,  this  routine  is  also  a linear  interpolator. 
Program  modifications  can  be  made  to  increase  the  number  of 
points  and/or  the  degree  of  the  polynomial  used  in  the 
fitting  procedure s» 


4d 


Figxire  24.  Flow  chart  tor  signal  analysis  program. 
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GNFILT: 

This  routine  will  calculate  the  coefficients  A(J)  and  B(J) 
for  a Butterworth  linear  recursive  low-pass  digital  filter. 
If  Y (K)  and  U (K)  are  the  output  and  input  arrays  to  the 
filter,  then 

N N+1 

Y(K)  = - Y.  A(J)*Y(K-J)+  ^ B(J)*U(K-J+1)  . 

J=1  J=1 

FILTER: 

After  calculating  the  coefficients,  the  routine  filters  the 
data. 

This  routine  filters  the  data  through  a Butterworth  linear 
low-pass  digital  filter  when  the  coefficients  are  known. 

SRTFUR: 

This  routine  performs  a bit  reversal  of  the  original  Nyquist 

m 

array.  If  I is  the  element  of  the  array,  then  I = 

n=o 

a(n)*(2**n),  a(n)  = 0 or  1.  The  element  specified  by  I is 

m 

exchanged  with  the  element  specified  by  J,  where  J = 

n=o 

a (m-n)  (2**n) . 

FFT: 

This  routine  calculates  the  Fourier  transform  of  a 
Nyquist  sampled  array  that  has  its  elements  in  bit-reversed 
order.  The  technique  is  generally  known  as  the  "fast 
Fourier  transform^" 

FWAL: 

This  routine  calculates  Walsh  transforms  of  a Nyquist 
sampled  array ^ The  technique  is  generally  known  as  the 
"fast  Walsh  transform^" 

AUTCOR: 

This  routine  calculates  the  autocorrelation  function  for  a 
Nyquist  sampled  array.  From  the  autocorrelation  function, 
the  power  spectrum  is  calculated  by  the  method  of  refined 
spectral  densities. 

PLOTS: 

This  routine  plots  an  equispaced  array. 

SAVE: 

Utility  subroutine  for  saving  arrays. 

UNSAVE : 

Utility  subroutine  for  retrieving  saved  arrays. 
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3.2  Input  Data  Descriptions 


Listed 
CARD  1: 


CARD  2: 


CARD  3: 


CARD  4: 


below  are  the  input  cards  and  their  format. 

READ  MTRACE 
FORMAT  (15) 

MTRACE:  The  total  number  of  the  sets  of  data  to  be 

processed  by  the  program. 

READ  FMT 
FORMAT  (20A4) 

FMT:  An  array  that  contains  the  format  of  the  input  data 

arrays  (X(I),  Y(I)). 

READ  TITLE 

FORMAT  (20A4) 

TITLE:  An  array  that  contains  80  characters  of 

alpha-numeric  information  describing  the  data. 

READ  NROT,  NSTAR,  NPOW,  ND^  JD. 

FORMAT  (515) 

NROT:  If  0,  then  the  input  array  must  be  rotated. 

NSTAR:  The  number  of  points  in  the  equispaced  Nyquist 

array.  It  must  be  an  interger  power  of  two. 

NPOW:  NSTAR  = 2**  NPOW. 

ND:  Used  for  plotting  power  spectra,  every  ND^^  point  will 

be  plotted. 

t h 

JD:  Used  for  plotting  time  series,  every  JD  point  will  be 

plotted. 
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CARD  5:  READ  INTER? , IFFFT,  IFFWAL,  IFAUTO,  IFGOF 

FORMAT  (615) 


INTER? 

= 

0 

Main 

calls 

NYQST 

= 

1 

Main 

calls 

LNYQ 

IFFFT 

= 

0 

Main 

skips 

FFT 

= 

1 

Main 

calls 

FFT 

IFFWAL 

= 

0 

Main 

skips 

FWAL 

= 

1 

Main 

calls 

FWAL 

IFAUTO 

= 

0 

Main 

skips 

AUTCOR 

= 

1 

Main 

calls 

AUTCOR 

IFGOF 

= 

0 

Main 

calls 

FILTER 

— 

1 

Main 

calls 

GNFILT, 

CARD  6:  READ  YSCAL,  XSCAL 

FORMAT  (2(2X,  E13.7)) 

YSCAL:  A scale  factor  that  multiplies  the  input  amplitude 

array, 

XSCAL:  A scale  factor  that  multiplies  the  input  time 

values. 

Note:  The  program  assumes  the  time  array  is  in  nanoseconds 

after  it  is  scaled  <> 

CARD  7:  READ  JF 

FORMAT  (2X,  14) 


JF:  The  number  of  points  in  the  input  time  and  amplitude 

arrays. 


CARD  8: 


READ  (XF(I),  YF(I),  1=1,  JF) 


CARD  9: 


CARD  10: 


CARD  11: 


CARD  12: 


CARD  13: 


FORMAT  (FMT) 

XF:  An  array  which  contains  the  initial  time  array. 

YF:  An  array  which  contains  the  initial  amplitude  array. 
READ  XZ  YZ  FORMAT  (2  (IX,  Ell. 4)) 

XZ,YZ:  Are  the  tablet  coordinates  of  the  zero  point  of  the 

trace. 

READ  XI,  X2,  Yl,  Y2 
FORMAT  (4F510) 

XI,  X2:  Coordinate  labels  for  digitization  points  on  the  X 

axis  of  the  grid. 

Yl,  Y2:  Coordinate  labels  for  digitization  points  on  the 

Y axis  of  the  grid. 

READ  (XT(I),  YT(I),  1=1,4) 

FORMAT  (2 (IX,  Ell. 4)) 

XT,YT:  Tablet  coordinates  for  the  four  grid  measurements 

defined  by  XI,  X2,  Yl,  and  Y2. 

READ  T,  V 

FORMAT  (2  (IX,  Ell. 4)) 

T:  A time-scale  factor  nanoseconds/division. 

V:  A voltage  scale  factor  volts/division. 

Note:  Cards  9,  10,  11,  and  12  only  appear  if  NROT  =1. 

READ  FI,  F2,  GW2 
FORMAT  (3E13.7) 

FI:  Half-power  point  for  the  low-pass  digital  filter  given 
in  megahertz. 
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F2:  <F1  is  another  frequency  at  which  we  specify  the  gain, 

GW2,  and  thus  determine  the  number  of  poles  for  the 
filter.  F2  is  in  megahertz  and  0 < 0.5. 

Note:  Card  13  only  appears  if  IFG0F=1.  If  IFGOF=0  then 

CARD  13:  READ  NP 

FORMAT (IS) 

CARD  14:  READ(A(D),  K=l,  NP) 

FORMAT  (5(2X,  E13.7)) 

READ  (B(K),  K=l,  NP+l) 

FORMAT  (5(2X,E13.7) ) 

If  further  analysis  is  desired  repeat  cards  2 through  12.  The  program 
will  continue  analyzing  data  until  all  input  data  are  exhausted.  If 
much  similar  data  are  to  be  processed,  it  would  be  convenient  to  define 
all  the  parameters  internal  to  the  program  and  then  delete  a great  many 
input  cards. 

A computer  printout  of  the  signal  analysis  program  is  included  as 
appendix  A. 

4.  REPRESENTATION  OF  EMP  WAVEFORMS  BY  PARAMETRIZED  FUNCTIONS 

As  previously  reported,  it  is  very  convenient  to  characterize  an 
EMP-induced  waveform  by  a finite  set  of  parametrized  functions.  This 
allows  an  analyst  to  conveniently  handle  the  thousands  of  waveforms 
generated  in  any  given  experiment.  At  present,  most  analysts  deal  with 
the  data  as  a digital  record  consisting  of  '\^10^  digital  values  per 
waveform  and  an  equivalent  number  of  digital  values  in  the  Fourier  power 
spectrum.  Although  this  is  a valid  approach,  it  is  difficult  and  tedious 
to  make  comparisons  between  large  sets  of  data  and  to  discover  trends  in 
the  data. 

A large  fraction  of  the  data  that  is  generated  in  EMP  experiments 
consists  of  waveforms  that  can  be  described  as  a product  of  a growth 
factor,  an  exponential  damping  factor,  and  a sinusoidal  function.  In 
many  of  these  records  there  is  more  than  one  dominant  frequency.  Most 

data  are  then  fitted  to  a judicious  mix  of  the  functions  where 

= t^  ^ exp(-C  t)sinu)  t . (75) 

n ^ n n ' 
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A linear  combination  of  the  generally  suffices,  thus  a function  ^ is 

fit  to  the  digital  data,  where 


N 

4>  = V a <J> 
^ n n 
n=0 


(77) 


or 


\ 


$ = y a $ 
^ n n 


N* 

E 


n=o 


b , 
n * n’ 


(78) 


There  is  no  general  rule  that  allows  a blind  selection  of  such  functions; 
thus,  at  some  point  the  data  must  be  examined  and  a "guess"  made  for  a 
good  set  of  functions  to  characterize  the  data. 

Following  is  a description  of  the  method  used  to  fit  an  N parcimeter 
nonlinear  function  to  a set  of  data.  Also,  some  examples  of  data  are 
included,  which  were  processed  by  using  the  computer  codes  especially 
developed  for  this  problem.  A listing  of  the  computer  program  with 
detailed  annotations  is  given  in  appendix  B. 

4.1  Method  of  Nonlinear  Least  Squares 

For  simplicity,  consider  the  problem  of  fitting  an  N-parameter 
function  with  one  independent  variable  X^,  $(X.,Pj  • . .P  ) , tg^^the  measured 
quantities  Y (X . ) where  the  subscript  i refers  to  the  i data  point. 
Thus,  we  find  tfte  parameters  that  minimize 


where^^  is  the  number  of  data  points,  and  is  the  statistical  weight  of 
the  i data  point. 

Let  us  assume  now  that  an  iterative  procedugg  has  been  defined  fg^ 
determining  the  parameter  P.;  that  is,  P.  = the  k iteration  of  the  j 
parameter.  Expand  the  function  in  ^ Taylor  serie|^  about  P.  and 
truncate  all  but  the  linear  terms.  ^^e  the  k iteration  to 

determine  the  parameters  for  the  (k+1)  iteration.  Let 
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(80) 


where 


and 


3 3 3 


k k+1 

Note  that  given  P . , we  must  now  determine  P . . Substituting 

equation  (80)  into  equation  (79)  yields  in  the  linear  Approximation 


s2~s‘  = 


(81) 


To  minimize  , 
k 


form 


2 

^fk  = 0,  Jl=l,...,N, 


(82) 


which  yields 


M 

Z 

i=l 


N 

Z 

j=i 


ji=i 


,N.  (83) 
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Defining  » ~ 


’i  " bhi 


C^  = ^ Z„.  W.  ; 

£ .^.  1 Jll  1 

1=1 


(84) 


then  from  equations  (81)  and  (82), 
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\ 1 / \ X i/ 

(85) 


f 

Or,  more  concisely. 


^ • AP  = C o 


(86) 


Equation  (83)  can  now  be  inverted  to  solve  for  the  AP 


e AP . hence, 

, thj  P . should 


the  P • If  now  our  iterative  procedure  is  converging, 
be  closer  to  the  values  P’^  which  minimizes  and  can  then  be  used  for 
the  iteration.  This^ iterative  procedure  is  continued  unti^  the  use 

of  P.  produces  a chi-squared  which  differs  from  that  using  P.  by  less 
than^some  preset  value,  V usually  V = 10"^;  that  is,  ^ 


"k+1 


-Sjl< 


V 


(87) 
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To  insure  that  the  step  AP . is  converging,  first  note  that  in  the  linear 
approximation  the  quantity^ 


-k  it  ,k . k ^ 


N N 


tl  ‘ i-1  j=l  «-l  J 1 » 


(88) 


M 

V W. 
i= 1 ^ 


[S  ^5^ 


k k 
AP . 
D 


> 0 


which  implies  the  must  always  be  positive  semidef inite.  Thus,  if  at 
the  k^  iteration  D <0,  then  change  the  sign  of  all  the  AP^.  Another 
approach  that  shows  this  is  to  expand  in  a Taylor  series  and  to  keep 
only  the  first  deriviatives  of  that  is. 


= S? 


Ik 

3S^ 

k AT.k  ^ 1 
APrt  + 

1 2 

3 

“ 8P„3P. 
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«-D  ^ 3 
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AP  AP . + ...  , 
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= W. 

i=i  ^ 


(<-.  ) - 2*.  Z Z Zji  j 

Aj  D ^ 


(89) 


.k  . 


If  AP^'  is  considered  as  defining  a vector,  cgnstr^ct  to  be  a function 
of  one  variable  a,  for  example,  by  letting  AP^-^uAP^;  then  we  have 


s2  (a)  = s2  - 2ttD^  + 


(90) 


Now  the  slope  of  (a)  evaluated  at  a=0  is 


dS^ (g) 


da 


~ ~ 2D„ 

a=0  T 


(91) 
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and  to  reacn  a minunum  tnis  slope  must  be  negative,  that  is,  ^ 
in  equation  (86). 

Although  the  condition  insures  convergence,  barring  roundoff 

errors,  the  procedure  can  take  divergent  stex:>s,  thus  oscillating  and 

giving  slow  convergence.  The  following  test  has  been  used  to  overcome 

this  problem  and  found  to  work  satisfactorily  on  most  data.  It  is  first 
necessary  to  check  and,  if  it  is  negative,  to  change  the  sign  of  ail 
the  Then  test  to  determine  if  the  value  if  it  is, 

perform  the  test  described  by  equation  (85)  ^ If  ,[(ot)  = 1]>  , 

, (oi)  must  be  calculated  for  the  following  four  values  of  Thus, 


a = 1/2  , 


[determined  by  the  parabola  with  slope  -D  at  a=0  and  passing  through  the 
points  s2  and  (a=l)  ] , 


a = l/o!j 

(determined  from  the  parabola  given  in  equation  (90) ) , 

“ * + s2^^(c.1)-2s2^P<.-1/2)] 

[determined  by  the  parabola  passing  through  the  points  S^,  . (a=l/2) , 

and  S^^^(a=l)]  with  the  restriction  that  any  value  of  a < 10  ^ is 
ignored.  From  these  four  (or  less)  values  of  S (a) , find  the  minimum, 

K.  I X 

compare  it  with  S^,  and  perform  the  termination  test  equation  (85)  if 
is  improvedo  If,  however,  is  still  not  improved,  then  either  bad 

data  points  remain  or  the  initial  starting  parameters  are  too  far  from 
the  ones  that  minimize  S^. 
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After  achieving  a minimum  as  determined  by  equation  (87)  , an  error 
analysis  is  performed  to  give  an  indication  of  how  well  the  parameters 
can  be  determined  from  the  data«  The  method  used  is  that  described  by 
Cohen,  Crowe,  and  Dumond.^  Their  discussion  deals  only  with  linear  least 
squares  fitting  and  we  have  not  made  a detailed  study  for  the  nonlinear 
case*  However,  in  the  neighborhood  of  the  minimum  in  the 
linearization  given  by  equation  (81)  should  be  a good  approximation * 
With  the  above  warning  the  error  (o^)  on  the  parameter  is  given  by. 


where 

>•  =s2/(M-N), 


(92) 


(93) 


which  is  the  chi~squared  normalized  with  the  number  of  degrees  of 
freedom^  To  be  conservative  in  our  error  estimate  is  set  equal  to 
one  if  the  fit  gives  a smaller  value. 

The  other  important  quantities  in  the  error  analysis  are  the 
correlation  coefficients  defined  by 


^3 


= (A-l) 


3j 


(94) 


The  correlation  coefficients  are  necessary  for  computing  effects  of 
error  propagation  when  using  the  "best**  parameters.  Consider  a function 
f which  depends  on  L fitted  parameters,  then  the  error  on  f is  given  by 


^f 


y 2f 

2P. 


2r 


(95) 


Ru  Conen,  Ku  Crowe,  and  J.  W.  Dumond,  Fundamental  Constants 
of  Physics,  Interscience  Publishers,  Inc.,  New  York  (1957),  Ch.  7. 
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4o2  Parametrization  of  an  EMP  Waveform  and  its  Autocorrelation  Function 


As  simple  examples  of  the  use  of  the  program  the  data  of  figure  4 
has  been  analyzed  in  several  different  ways.  The  power  spectrum  of  the 
time  series  data  is  shown  in  figure  5.  It  is  seen  that  there  are  two 
predominant  frequencies  in  the  time  waveform^  13.2  and  24  MHz.  To 
isolate  one  of  the  frequencies  (and  obtain  a less  complicated  time 
series) , the  data  of  figure  4 were  digitally  filtered,  using  a B-pole 
low-pass  digital  filter  with  a half-power  point  at  17.5  MHz  and  a gain 
at  25  MHz  of  0.05.  The  resulting  time  series  is  shown  in  figure  25, 
Its  power  spectrum  is  shown  in  figure  26.  It  is  seen  from  figures  25- 
and  26  that  the  time  series  consists  primarily  of  a sinusoid  with 
frequency  of  13.2  MHz.  To  the  time  series  of  figure  25  is  fit  the 
following  function 


^(t,P)  = P(l)  t exp  (-P(2)  t)sin[P(3)  t P(4)]  (96) 

It  took  the  program  five  iterations  to  converge  to  a solution.  The 
initial  estimates  and  the  final  fitted  parameters  are  given  in  table  II. 
A combined  plot  of  the  experimental  data  and  the  fitted  function  is 
shown  in  figure  27o  Generally,  the  fit  is  good  only  in  the  central 
portion  of  the  trace  and  it  fails  badly  at  the  beginning  and  end  of  the 
traceo  This  clearly  means  that  equation  (96)  is  not  the  best 
representation  of  this  trace « 

From  the  filtered  time  series  of  figure  25,  the  autocorrelation 
(lagged  products)  function  was  formed  as  displayed  in  figure  28.  The 
following  function  was  fit  to  these  data. 

4>(t,P)  = P(l)exp(-P(2)  t]cos[P(3)  t]  (97) 


It  took  the  program  two  iterations  to  converge  to  a solution.  The 
initial  estimates  and  the  final  fitted  parameters  are  given  in 
table  III.  A combined  plot  of  the  data  and  the  fitted  functions  is 
shown  in  figure  29.  As  can  be  seen  the  fit  is  quite  good  over  the 
entire  trace. 

As  our  third  example,  the  autocorrelation  function  for  the  data  of 
figure  4 was  calculated  and  is  plotted  in  figure  30  (no  filtering  of  any 
kind  was  done) ; it  was  fit  by  the  following  function. 


4>(t,P)  P (1)  exp(-P  (2)  t]  cos  [P  ( 3)  t]  + P (4)  exp  [-P  (5)  t]  cos  [P  (6)  t] 


(98) 
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Figure  25.  Time  series  of  figure  4 after  digital  filtering. 
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Figure  26.  Power  spectrum  for  the  time  series  of  figure  6. 

It  took  the  program  three  iterations  to  converge  to  a solution. 
The  initial  estimate  and  the  final  fitted  parameters  are  given  in 
table  IV.  A conbined  plot  of  the  experimental  data  and  the  fitted 
function  is  shown  in  figure  31.  The  fit  is  cjuite  excellent  throughout 
the  length  of  the  trace. 

Several  attempts  were  made  to  fit  the  data  of  figure  4 with 
various  forms  of  the  functions  listed  in  equation  (76) . Whenever  n was 
greater  than  2,  the  fitting  procedure  either  diverged  or  produced  a bad 
fit. 
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TABLE  II.  INITIAL  ESTIMATES  AND  FINAL  FITTED  PARAMETERS  FOR  EQUATION  (96) 
AND  DATA  OF  FIGURE  25 


Parameter 

Initial  estimates 

Final  fitted  values 

p(i) 

1.08 

+0.742 

P(2) 

4.27 

+4.95B 

P(3) 

90,43 

+84.433 

P(4) 

-19.59 

-18.207 

TABLE  III.  INITIAL  ESTIMATES  AND  FINAL  FITTED  PARAMETERS  FOR  EQUATION  (97) 
AND  DATA  OF  FIGUPJ?  28 


Final  fitted  values 
+0.00102 
+3.921 
+84.196 

TABLE  IV.  INITIAL  ESTIMATES  AND  FINAL  FITTED  PARAMETERS  FOR  EQUATION  (98) 
AND  DATA  OF  FIGURE  30 


Parameter  Initial  estimates  Final  fitted  values 


p(i) 

0.0014 

O0OOIO8 

P(2) 

3.92 

4.0289 

P(3) 

00 

c 

to 

o 

84.7402 

P(4) 

0.001 

0.00148 

P(5) 

3.92 

9.8083 

P(6) 

147.58 

146.0411 

Parameter 

Initial  estimates 

P(l) 

+0.001 

P(2) 

+4.27 

P(3) 

+86.35 
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AMPUTUOe 


Figure  27*  Plot  of  a least-squares  fit  of  equation  (96)  to  the  time 
series  of  figure  6. 
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Figure  28.  Autocorrelation  function  for  the  time  series  of  figure  6. 
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Figure  29.  Plot  of  a least-squares  fit  of  equation  (97)  to  the  auto- 
correlation function  of  figure  21. 
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Figure  30.  Autocorrelation  function  for  the  time  series  of  figure  1. 
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It  seems  clear  that  before  fitting  a function  to  the  data  (to 
obtain  parametrization)  a process  such  as  digital  filtering  or  formation 
of  the  autocorrelation  function  is  necessary  to  isolate  the  functional 
dependence  (at  least  the  oscillatory  part) - 

Following  are  parameters  to  the  power  spectrum  of  the  record. 
Since  the  power  spectrum^  P(cj)  , is  given  by  the  cosine  transform  of  the 
autocorrelation  function,  we  have 

CO 

P(U))  = 2 I cosa)t[P(l)exp{-P(2)  t}cos[P(3)  t] 

•f  P(4)exp[-P(5)  t]cos{P(6)  t}]dt  . (^^3 

Thus,  after  some  algebraic  manipulation,  a power  spectrum  is  given  by 


^ P(l)  P(2)  ^ P(4)  P(5)  ^ (100) 

[o)-P(3)  ]2-fP(2)2  [o)-P(6)  ]2+P(5)2 


It  would  have  been  just  as  easy  to  fit  equation  (100)  to  the  power 
spectrum  displayed  in  figure  5 and  thus  obtain  a parametrization  of  the 
data.  However,  we  have  used  a plot  of  equation  (100)  superimposed  on 
the  power  spectrum  of  figure  5 as  shown  in  figure  32. 

Considering  the  storage  of  data,  it  is  concluded  that  only  six 
numbers  must  be  stored  to  characterize  the  pulse  instead  of  thousands. 
Thus,  long-term  storage  costs  are  minimized  and  a high  degree  of  \ 
analytical  ease  is  realized  in  handling  large  blocks  of  data. 

4.3  Program  Listing  and  Description 

The  program  is  written  in  FORTRAN  and  consists  of  eleven 
suiiroutines  and  two  user  supplied  functions;  the  program  names  with  a 
description  of  each  are  listed  below. 

MAIN  This  routine  directly  or  indirectly  calls  the  rest  of  the 
subroutines  and  thus  controls  the  passage  of  the  program 
through  all  the  calculations  used  in  the  fit. 

READIN  This  routine  reads  in  all  the  input  parameters  and  data  as 
well  as  setting  up  all  the  internal  control  parameters. 

LSQPHI  This  routine  calculates  the  theoretical  function,  residuals, 
and  chi  square.  The  user  must  supply  the  function  PHIFNC(I) 
for  LSQPHI  to  use. 
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o PLOT  OF  THE  POWER  vs 
FREQUENCY  CALCULATED 
BY  THE  METHOD  OF  RE- 
FINED SPECTRAL  DENSITY 

— PLOT  OF  THE  POWER  vs 
FREQUENCY  CALCULATED 
FROM  EON  100 
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FREQUENCY  (UNIT  0.5989  MHz) 


Figure  32.  Plot  of  power  spectrum  obtained  from  equation  (100)  and  that 
obtained  from  a numerical  algorithm. 


ELTS(I)  This  routine  calculates  all  the  derivatives  of  the 

theoretical  function  at  the  data  point.  The  user  must 
supply  the  function  (DPHI(I)J)  for  ELTS  to  call. 

GRIND  This  routine  sets  up  the  matrix  of  derivatives,  the 

constants  in  the  normal  equations,  then  calls  the  matrix 
inverter. 
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TESTS 

This  routine  does  the  necessary  testing  of  the  changes  in 
tlie  parameters  and  performs  the  necessary  changes. 

ERRMAT 

This  routine  does  all  the  calculations  for  a full  error 

analysis  on  the  fitted  parameters. 

CALS2  This  routine  is  called  by  TESTS  repeatedly  for  the  purpose 


of  swapping  arrays  and  computing  chi. 

INVMAT 

Inverts  a symmetric  N x N matrix. 

SCRIBE 

This  routine  writes  out  all  the  results  of  the  program  fit. 

PLOTS 

Gives  a plot  of  the  experimental  data  and  theoretical 
function  on  the  same  graph  so  that  the  goodness  of  the  fit 
can  be  clearly  seen 

PHIFMC(I) 

PHIFNC  = 0(X(I),  P(l),  P (2) , . . . ,P (N) ) 

DPHI  (I ,J) 

DPHI  = tIttT  (X(I),  P(1),  P(2)  , . . . ,P(N)  ) 
oP  ( J) 

The  flow  chart  of  figure  33  shows  the  logical  connection  of  the 
subroutines  and  functions. 

Listed  below  are  the  input  parameters,  their  meaning,  and  their 


format. 

CARD  1 

READ  JSTOP,FMF  FORMAT  (I5,5A4) 

JSTOP  Maximum  number  of  iterations  allowed  for  each  fit. 

This  is  used  to  prevent  excessive  iterations  in  case 
there  are  errors  in  the  experimental  data  or  control 
parameters.  Usually  JSTOP  100. 

FMT  The  format  with  which  the  experimental  data  are  to 
be  read  in,  e.g.,  (6E13.7). 

NOTE:  Card  1 appears  once  independent  of  the  number  of 

sets  of  data  to  be  fit. 

CARD  2 

READ  TITLE  FOR^IAT  (20A4) 

TITLE  Any  combination  of  alphanumeric  data  may  appear  in 
columns  1 through  80  to  identify  the  fit. 
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Figure  33.  Flow  chart  of  nonlinear  least-squares  fitting  program. 
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CARD  3 READ 

FORMAT 

NCHAN,  NBADPT,  NRUNS,  NRP,  NPSTEP,  IFPLOT,  IFCOR 
(715) 

NCHAN 

The  nuinber  of  experimental  data  points. 

NBADPT 

The  number  of  bad  points  in  the  data. 

NRUNS 

The  number  of  sets  of  experimental  data  to  be  added 
together. 

NRP 

The  number  of  independent  parameters. 

NPSTEP 

^0,  used  for  plotting;  that  is,  every  NPSTEP^^  data 
point  is  plotted  starting  with  the  first  (0  is 
internally  set  to  1) . 

IFPLOT 

= 0 no  plot  of  data  and  fit 
= 1 gives  plot 

IFCOR 

= -1  This  prints  and  punches  the  correlation 

coefficients 

= 0 Does  not  print  or  punch 
= 1 Only  prints 

CARD  4 READ 

(P(T) , 1=1,  NP) 

P 

the  array  containing  the  fitted  parameters. 

CARDS  5 READ 

(Y(I),  1=1,  NCHAN)  FORMAT  (FMT) 

Y 

the  array  containing  the  experimental  data,  one  run 
following  another  for  NRUNS  worth  of  data. 

CARDS  6 READ 

(X(I),  1=1,  NCHAN)  FORMAT  (FMT) 

X 

the  array  containing  values  of  the  dependent 

variable  for  ^ and  Y;  for  example,  0(X(I),P). 

CARDS  7 READ 

(W(I),  1=1,  NCH7VN)  FORMAT  (FMT) 

W 

the  array  conta^j^ing  values  of  the  statistical 
weights  for  the  I data  point. 
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CARDS  8 READ 


(NBAD(I),  I=Jr  NBADPT) 


NBAD  Index  number  of  the  bad  data  points  in  any  order. 

These  cards  do  not  appear  in  the  data  deck  if 
NBADPT  = 0. 

If  multiple  fits  are  desired,  repeat  cards  2 through  8.  The 
program  will  continue  doing  fits  until  it  has  exhausted  all  of  the  input 
datao 
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APPENDIX  A.  SIGNAL  ANALYSIS  PROGRAM 


( JAN  73  ) 


OS/360  FORTRAN  H 
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COMPILER  OPTIONS  - NAHE=  MA I N ,OP T*02 # LI NECNT*57 , S I ZE*3D03K, 

SOURCE fEBCOICtNOL 1ST tOECKf LOAD* MAPtNOEOITt IO»NOXREF 

C 

MAIN  PROGRAM  FOR  THE  ANALYSIS  OF  DIGITAL  TIME  SERIES 
PROGRAMMER  OR.  THOMAS  A TUMOLILLO 

USAMC  HARRY  DIAMOND  LABORATORIES 
WASHINGT0N80.C.  20A38 


C3MMON  VARIABLES  DEFINITION 


YNYQA 


NW 

NR 


NT 

XF 

Y«= 

NSTAR 

NPOW 

NO 

JF 

CT 

ST 

XO 

YO 

xs 

YS 

T 

V 

XNYO 


COMPLEX  ARRAY  USED  TO  STORE 
STAGES  OF  THE  ANALYSIS. 


THE  DIGITAL  TIME  SERIES  AT  VARIOUS 


SYMBOLIC  DESIGNATOR  FOR  THE 
SYMBOLIC  DESIGNATOR  FOR  THE 
SYMBOLIC  DESIGNATOR  FOR  THE 
ARRAY  WHICH  INITIALLY  CONTAINS 
ARRAY  WHICH  INITIALLY  CONTAINS 


OUTPUT  LINE  PRINTER 
INPUT  CARD  READER 
OUTPUT  CARO  PUNCH 
THE  TIME  VALUES 
THE  AMPLITUDE  VALUES 


INTERPOLATED  POINTS9IT  MUST  BE  A POWER  OF  TWO 


THE  NUMBER  OF 
NSTAR=2**NPCW 

USEO  FOR  PLOTTING  POWER  SPECTRA3EVERY  NO-TH  POINT  WILL 
THE  NUMBER  OF  POINTS  IN  THE  INPUT  DIGITAL  TIME  SERIES 
THE  COSINE  OF  THE  ROTATION  ANGLE 
THE  SINE  OF  THE  ROTATION  ANGLE 
X-COORDINATE’  OF  THE  GRATICULE  ORIGIN 
OF  THE  GRATICULE  ORGIN 

TABLET  INTEGERS/SCALE  DIVISION  ALONG  X AXIS 
TABLET  INTEGERS/SCALE  DIVISION  ALONG  Y AXIS 
NANOSECONDS/SCALE  DIVISION 
VOLTS/SCALE  DIVISION 

SAMPLING  INTERVAL  OR  THE  TIME  INTERVAL  BETWEEN 


BE  PLOTTED 


Y-COORDINATE 
SCALE  factor 
SCALE  FACTOR 
SCALE  FACTOR 
SCALE  FACTOR 
THE  NYOUI ST 


oooooooonnooooo 


f I 


BETWEEN  DIGITAL  VALUES 

XZ  THE  X-COORO  OF  THE  ZERO  POINT  OF  THE  TRACE 

YZ  THE  Y-COORO  OF  THE  ZERO  POINT  OF  THE  TRACE 

JO  USED  FOR  PLOTTING  TIME  SERIEStEVERY  JO-TH  POINT  WILL  BE  PLOTTED 

INTERP  INTERP=0  CALL  NYQST,=l  CALL  LNYQ 
IFFFT  IFFFT»0  NO  FFT,=l  CALL  FFT 
IFFWAL  IFFWAL-0  NO  WALSH  TRANSFORMt=l  CALL  FWAL 
I=AUT3  IFAUTO=0  NO  FT  VIA  AUTOCORRELAT I ONt = I CALL  AUTCOR 
IFGOF  IFG0F=0  CALL  FILTER,*!  CALL  GNFILT 

IFPUNF  IFPUNF*!  PUNCH  OUT  FREQUENCY  AND  POWER  SPECTRA  VALUES, =0  DONT 
IFPUNW  IFPUNW*!  PUNCH  OUT  SEQUENCY  AND  WALSH  POWER  VALUES, *0  DONT 
FPUNA  IFPUNA*!  PUNCH  OUT  TIME  AND  AUTOCORRELATION  VALUES, *0  DONT 

IFPUND  IFPUND*!  PUNCH  OUT  TIME  AND  AMPLITUDE  VALUES  OF  INPUT  DATA,=0  DON 

If=PUNI  IFPUNI  = l PUNCH  OUT  TIME  AND  INTERPOLATED  VALUES, *0  DONT 
IFPUNR  IFPUNR=l  PUNCH  OUT  FREQUENCY  AND  REFINED  SPECTRAL  DENS  I T I ES,*ODON 
COMMON  YNYQA(2O^0» 

COMMON  NW,NR,XF (2048) , YF ( 2046) ,NST AR ,NPOW ,ND , JF , CT , S T , XO, YO, 

I XS,YS,T,V,XNYQ 
2,NT,XZ,YZ, JD 

COMMON  INTERP, IFFFT, IFFWAL, IF AUTO, IFGOF, IFPUNF, IFPUNW, I FPUNA, 
IIFPUND,IFPUNI , IFPUNR 
COMPLEX  YNYQA 
NW  = 6 
NR*5 
NT*7 

READ  (NR,l)MTRACE 
I F3RMATII5) 

DO  120  I*l,MTRACE 
CALL  REAOIN 
CALL  CSTOUT 
IFIINTERP)  10,10,20 
10  CALL  NYQST 

GO  TO  30 
20  CALL  LNYQ 


■i 


03 


30  CALL  PLOTSUI 
CALL  SAVE(0» 

NPOWS=NPOW 
NK6EP=NSrAR 
IFIIFGOFJ  40,40,50 
40  CALL  FILTER 

CALL  SAVE(l» 

G3  ro  60 

50  CALL  GNFILT 

CALL  SAVF.m 
60  CALL  PL0TSI4) 

IFdFFFTI  80,80,70 
70  CALL  SRTFUR 

CALL  FFT 
CALL  PL0TSI2I 

80  IFdFAUTOI  100,100,90 

90  CALL  AUTCORd  ,NKEEP,NPCWSI 

100  IFdFFWAL)  120, 120, UO 

dO  CALL  UNSAVEIU 

CALL  FWAL 
CALL  PL0TS(2I 
120  CONTINUE 

STOP 
END 
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COMPILER  OPTIONS  - NAME=  MA I N, OPT=02 »LI NECNT=57 , S I ZE=03D0K » 

SOURCEtEBCOICtNOL 1ST « DECK* LOADt HAP* NOEDITt IO*NOXREF 
SUBROUTINE  REAOIN 
COMMON  YNYQAI20A8I 

COMMON  NH*NR,XF(20A8»  *YF(2O40I ,NSTAR*NPOW,ND»  JF, CT , ST , XO* YO , 

I  XS* YStTtVfXNYO 
2,NT,XZ,YZt JO 

COMMON  INTERP*IFFFT,IFFWAL,IFAUTO,IFGOF, IFPUNF, IFPUNW* IFPUNA* 
lIFPUND,IFPUNIt IFPUNR 
COMPLEX  YNYQA 

DIMENSION  FMTI20I ,TITLE(20I 
READINR* IIFMT 

1 F0RMATI20AAI 
READINRrl)  TITLE 
WRITE(NW,2)  TITLE 

2 FORMATllHl ,20X,20A4) 

READINR, 3 » NROT,NSTAR ,NPOW ,ND, JD 

3 F0RMATI5I5I 
W^ITE  (NW,13I 

13  FORMAT  llX,'NROT  NSTAR  NPOW  ND  JD* I 
WRI TE (NW* 14)  NROT  *NSTAR*NPOW*ND*  JD 

14  F0RMAT(I4,4I5,/) 

RE  AD (NR, 3)  I NTERP , I FFF T , IFF WAL , I FAUTO, I FGOF 
WRITE(NW,ll) 

11  FORMAT! IX, 'INTERP  IFFFT  IFFWAL  I FAUTO  IFGOF') 

KRI TE (NW, 12) INTERP, IFFFT , IFFWAL , I F AUTO, I FGOF 

12  F0RMAT11X,2I6,2I7, I6,6I7,Z) 

READINR, 8)  YSCAL,XSCAL 

8 F0RHAT(2I2X,E13.7)  ) 

WRITE  (NW,15) 

15  FORMATIIX,*  Y SCALE  FACTOR  X SCALE  FACTOR  •) 

WRITE  INW,16)  YSCAL,XSCAL 

16  F0RMATI2(2X,El 3.7),/) 

WRITE(NW,5) 


5 FOPHAT(»  IMAGE  OF  INPUT  DATA’) 
RFADINR,7)  JF 

7 F3RMATI2X, 14) 

PE  ADINR ,FHT ) I XF I I ) , YF I I ) » I » JF ) 
no  100  1=1, JF 

YF I 1 ) = YSCAL*YF ( I ) 

100  XF (I )=XSCAL*XF( I) 

X'4Y0  = XF(  JF)/NSTAR 

WR1TE(NK,4)  I XFI  I ) , YF I I ) , 1 = 1, JF) 

4 FORMAT I9I IX,EI I .4 ) ) 

IF INROT-l )60, 50,50 
50  CALL  SCARTP 

CALL  ROT 
WR  ITE(NW,6) 

6 FORMAT!/, • SCALED  AND  ROTATED  DATA’) 

WRITE(NW,4)  IXF( I ) ,YF(  I ) , 1 = 1, JF ) 

50  RETURN 

END 
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COMPILER  OPTIONS  - NAME*  M A I N, 0PT=02 , L I NECNT  = 5 7 , S I ZE*0000K , 

SOURCE » EBCDIC fNOL 1ST ,OECK,LOAO» MAP, NOEO I T » lOf NOXREF 
SUBROUTINE  GNFILT 
COMMON  YNYQA(2048) 

COMMON  NH,NR ,XF(2348I , YF ( 2348 1 , NST AR, NPOW, NO, JF ,CT , ST , XO, YO, 

I XS,  YS,  T,  V,XNYO 
2,NT,XZ, YZ, JO 

COMMON  INTERP,IFFFT, IFFWAL,IFAUTO, IFGOF, IFPUNF, IFPUNH, I FPUNA, 
lIFPUNO,IFPUNl, IFPUNK 
COMPLEX  YNY0A,U,R ,S,ROOT,BETA,UNIT 
EQUIVALENCE  ( ARR AY ( I t , YNYQA ( LH 
DIMENSION  ARRAYI4096I 
DIMENSION  STHI43) ,CTH(43» , IN0M(43) 

DIMENSION  U(40) ,R(20) ,S(20) ,ROOT(20) ,P( 20) , 
lOI 20) ,PQ(20) ,A( 20) ,B( 20) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCI 

THIS  ROUTINE  WILL  CALCULATE  THE  COEFFICIENTS  AU)  AND  B(J) 

FOR  A BUTTERWORTH  LINEAR  RECURSIVE  FILTER.  IF  Y(K)  AND  U(K) 

ARE  THE  OUTPUT  AND  INPUT  ARRAYS  FOR  THE  FILTER  THEN, 

YIK)=-SUM(1  TO  N) (A( J)«Y(K-J) )^SUM(l  TO  Nf I ) I B (J ) «U ( K- J ♦ I ) ) 

THE  CAIN  FACTOR  FOR  THE  BUTTERWORTH  FILTER  IS  GIVEN  BY 

GAIN* /HI  W)  /♦♦2=I/(H-(TAN(W*T/2)/TAN(WI*T/2)  )♦♦(  2*N)  ) 

HERE  T IS  THE  SAMPLING  INTERVAL,  N IS  COMPUTED  BY  SPECIFYING 
THE  HALF  POWER  POINT  Wl  AND  THE  RELATIVE  GAIN  AT  THE  POINT 
W=W2>WI,  AND  THE  TRANSFER  FUNCTION  HIZ)  IS  GIVEN  BY  THE  TWO  FORMULAS 

H(Z)*SUMII  TO  NH)  (B(J)*Z**(-J*^l)  )/(  USUMI  I TO  N ) I A ( J ) ♦ Z**  I - J ) ) 

H(Z)=BETA*( UZ)**N/(Z-P( I ) ) ♦( Z-P 1 2 ) ♦. . •♦ ( Z-P IN ) ) . 
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c 

C HERE  THE  POLES  OF  THE  FILTER,  PII),  ARE  GIVEN  BY 
C 

C P(  I l = ( l-TAN(ARG)**2+l-l  I**95*2*rAN(ARGI*SIN(TH( 

C ( l-2*TAN(  ARGI*COSI  TH(  I )l  ♦•TANI  ARGI**2» 

C WHERE  ARG=Wl*T/2,TH(  n=n-ll*PI/N  FOR  N ODD,  = ( 2*  I- II  *P  I /2*N 
C FDR  N EVEN,  AND  ONLY  THOSE  N Ptll  ARE  CHOSEN  THAT  SATISFY 
C /PIIKl..  THE  NORMALIZING  FACTOR  BETA  IS  GIVEN  BY 

C 

C BE  TA=  ( l-P(ll  !♦(  l-PI  21  )$...♦(  l-P(NM/2**N 
C 

C BY  COMPARISON  OF  THE  TWO  EXPRESSIONS  FOR  THE  TRANSFER  FUNCTION 
C THE  FILTER  COEFFICIENTS  AUI  AND  B(J)  ARE  DEDUCED 
C 

CCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
PI =3. 1415927 

c 

C READ  IN  OR  SOMPUTE  THE  FREQUENCIES  FI  AND  F2,  THE  GAIN  AT  F2,GW28  AND 
C THE  SAMPLING  INTERVAL 
C 

TI ME=0,00l*XNYQ 
REAO(NR,l)  Fl,F2,GW2 

1 FORMAT  13613.7) 

WRITE(NW,2)  Fl,F2,GW2,TIME 

2 FORMATIlHl,'  THE  HALF  POWER  POINT  IS  • , E13.7, ‘MHZ' »/ , 

I*  THE  GAIN  AT  ',E13.7,'MHZ  IS*,E13.7,/, 

2*  THE  SAMPLING  INTERVAL  IS  *, E 1 3. 7, • MICROSECONDS ',/ I 

C 

C COMPUTE  AND  PRINT  OUT  THE  NUMBER  OF  POLES  FOR  THIS  FILTER 
C 

NMAX=lO 

ARGl»Fl*TIME*PI 
TAGl=TAN( ARGl I 


I 


I V 


A^G2=F2* TIME^P I 
? T = ( TAN( ARG2  J /TAG  1 1 **2 
TVAL=( I .0-GW2) /GW2 
N=  I 

N3E=-l 

13  IF  (RT«*,N.GE.  TVALI  GO  TO  20 

N=NU 
NQF=-NOE 

IF  tN.GT.NMAX>  GO  TO  25 
GO  TO  13 

25  WRITE(NW,3)  NMAX 

3 FORMATtlX,'  THE  NUMBER  OF  POLES  FOR  THE  FILTER  EXCEEDED  ',/♦ 

l»  THE  ALLOWED  LIMIT  OF • , I 5 I 
STOP 
C 

C COMPUTE  ALL  THE  2*N  POLES  WHICH  ARE  EITHER  INSIDE  OR  OUTSIDE  THE 
C UNIT  CIRCLE 
C 

23  TW0N=2*N 

I TWON=  TWON 
FN  = N 

WR I TE (NW,AI 
WR I TE (NW,5> 

03  ^0  I=l,ITWCN 
IFtNCE.EO.il  GO  TO  30 
TH=t I-I I *PI/EN 
G3  T3  35 

30  TH=t 2*1-1 >*PI/TWON 

35  CTHlII=COSlTHI 

STHt  I l = SINt  THI 
40  CONTINUE 

Rl=l-TAGl**2 
R2=2*TAGl 
01=1 ♦TAGl**2 


V xicmidv 


K=0 

DO  50  l=l,ITWCN 
DEN0M=Dl-R2*CTH(  I ) 

P( I >=R l/DENOM 

Q{ I )=R2*STH( I » /DENOM 

Pan  i = Pt  I l♦*2♦o(  I 

WRITE(NW,6)  I »P(1  ) rOn  i tPQ(  I i 

c 

C SELECT  ONLY  THOSE  POLES  WHICH  ARE  INSIDE  THE  UNIT  CIRCLE 
C 

IF  tPQI I » .GT. I.O)  GO  TO  50 
K = K«-l 

ROOT IKI=CMPLX( PI  I 1,01111 
50  CONTINUE 

5 FOKMATIIH/, • INDEX* ,5X,  ‘REAL  P AR I • , 5X , * I M AG  P ART • , 8X , • SOU AR E • I 

^ FORMAT! IH/,20X, *F ILTER  POLE  LOCATIONS'1 

6 F0RMAT(IX,I5,3(IX,EI3.7H 
IF (K.EQ.NI  GO  TO  60 
WRITEINW,?)  KtITWON 

7 FORMATIIH/,'  ONLY', 13,'  POLES  OUT  OF', 13,'  ARE  INSIDE  UNIT  CIRCLE' 
I,/, IX,'  PROGRAM  EXECUTION  HALTEO'I 

STOP 

C 

C COMPUTE  BETA 
C 

60  BETA=CMPLXll. 0,0.01 

UN  I T=BETA 
00  70  1=1, N 

70  BETA=BETA*(UNIT-ROOT( m /2.0 
C 

C COMPUTE  THE  BINOMIAL  COEFFICIENTS 
C 

I NOM I II  = l 
NP  I=N<-  I 
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CD 


D1  75  1=2, NPl 

75  INJOMd  > = ( (N-d2)*INOM(  I-l)  >/(I-U 

WRITE(NW,14»  N 

14  PnRMArC  BINOMIAL  COEFFICIENTS  FOR  I I Z I , I 2 , / , 
l»  INDEX  COEFFICIENT'! 

WRI TEINW, 15»  ( I , INOMI I ) , I=l ,NPl) 

15  FORMAT (4X, 12, 7X, 14! 

C 

C COMPUTE  THE  DIGITAL  FILTER  COEFFICIENTS 
C 

IF (N-2 180,90, 100 
80  Uin=-ROGTm 

GO  TO  125 

90  01 3I=:MPLX( 1.0,0.01 

U( 2)=-P30TI 1)-R00T(2I 
01  ll=R03Tm*RO0TI2» 

CO  TO  125 

100  Ul ll=R03TIll*R00TI2l 

U(  2)=-ROOTm-ROOT(2l 
Ul 3»=CMPL<( 1.0,0. 0) 

L = 3 

R12I=CMPLXI 1.0,0.01 
on  120  K=3,N 
PI 1»=-RC0TIKI 
33  I 10  I =l,L 
110  S(  n=0(  I I 

CALL  P0LMLTI2,R,L,S,LR,U» 

L = LR 

120  CONTINUE 

125  WRITE(NW,8» 

8 F3RMATIIH/,'  I NDE X ' , 5 X , • RE AL  A ( II • , 5X , • I MAG  AIII'I 

, 00  130  I=l,N 

IS=Ni-l-I 
RA=REALIU( ISI I 
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A(  I »=RA 

AA  = AIMAG(U(  ISn 
130  WRITE(NW,9)  I,RA,AA 
9 F0RMAT(IX,I6,2( lX,El3.7n 

NP l*Nf I 
WRITE(NM»m 

II  FORMATUH/,*  INDEX', 5X. 'REAL  d ( II  ' , 5X  , ' I MAG  8(1)') 
NP2=Nf2 

DO  lAO  1=1, NPl 

NR=NP2-I 

SHI  T = I NOM( NR) 

SI  I )=CMPLX( SHI T,0.0) *BETA 
R3  = REALI SI  I ) ) 

B( I )=RB 
33=AIMA3(S( I ) ) 
lAO  WRITF(NW,9)  I,RB,BB 

DIGITALLY  FILTER  THE  INPUT  TIME  SERIES 

DO  170  K=l,NSTAR 
Z = 0.0 

in  150  J=l,NPl 
nD=K-jf  I 

IFIIND.LE.O)  GO  TO  155 
2=Z»BI J) *AIMAG( YNYQAI IND) ) 

150  CONTINUE 
155  DO  160  J=l,N 
IND=K-J 

IF  I IND.LE.O)  GO  TO  165 
Z=Z-A( J) *REAL ( TNYOAI IMD) ) 

160  CONTINUE 

165  X=AIMAG( YNYUA(K) ) 

YNYQA(K)=CMPLX(Z,X» 

170  CONTINUE 


o o n 


DJ  180  K=l,NSIAR 
X=REAL(VNYGA(K) > 

180  YWYOAUI  =CyPLX(  X,0.0» 

NS2=2*NSTAR 

PRINT  OUT  AND  PLOT  THE  FILTERED  TIME  SERIES 
WO  I TE (NW,  12) 

WR I TE (NW, 131  I ARRAY  I N , I = l,NS2t 2» 

12  FIRMATC  FILTERFI)  TIME  SERIES*! 

13  FORMAT (3 ( I X, El  3. Til 
RETURN 

END 


VC 

c 
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COMPILER  OPTIOMS  - NAME=  MA I N , OP T =32 , LI NECNT =5 7 , S I Z F= DD33K , 

SOURCE ,EBCDIC,NOL  1ST  ,OECK, LOAD, map, NOEDIT, ID,NOXRFF 
SUBROUTINE  FILTER 

THIS  IS  A RECURSIVE  FILTER  ROUTINE  WHERE  NP  IS  THE  NUMBER  OF  POLES,  A 
AND  B ARE  THE  FILTER  COEFFICIENTS, 

AIMAGl YNVJA)  IS  THE  INPUT  ARRAY, 

REALIYNYOAI  IS  THE  OUTPUT  ARRAY 

RE  AL  ( YNYOA  ( I » )=SUMB(J)  *A  I M AG  ( Y NY  QA  (I  -JH  I - SUMA  ( J I *REAL  I YNYU  A I I- J N . 
FOR  NO  FILTERING  SET  B ( U =l  , B ( K » =3, K. NE . I ; A I K I = 3 ALL  K. 

COMMON  YNYQA(20A8» 

C3MM0N  NW,NR,XFI2  048»  , YF ( 2048  I , NS  I AR , NPOW , ND , JF , CT , ST , XO , YO , 

I XS,  YS,  T,  V,  XNY3 
2,NT,XZ,YZ, JD 

COMMON  INTERP,IFFFT, I FF W AL , I F AUTO , 1 FGOF , IFPUNF, IFPUNW, IFPUNA, 

I IFPUND , IFPUNl , IFPUNR 
COMPLEX  YNYQA 
DIMENSION  ARRAY(4396) 

DIMENSION  A(  ICI  ,B(10I 
EQUIVALENCE  I ARRAY!  I I , YNYQA  lU  ) 

READ(NR,3»  NP 

3 FORMATIISI 
READ(NR,4)  ( A ( K ) , K =l , NP I 
NP  l=NPf I 

RE  AI)(NR,4)  ( B!KI  ,K=l  ,NPl  ) 

4 FORMAT  151 2X,tl 3, 7) » 

WR I TE(NW,  lOl > 

101  FORMAT ( IHl , » FILTER  COEFFICIENTS'! 

WR ITEINW,  102!  ( K , A ( K I , K , B ( K ) , K= I , NP ) 

102  FORMAT  nx,'A(',I2,'!  = ',El3.7,'B(',I2,'I  = ',El3.7) 

WR1TE(NW,103»  NPl,B(NPll 

103  F0RMAT(21X,'B(',I2,')=’,E13.7) 

C DIGITALLY  FILTER  THE  INPUT  TIMESERIES 
C 


DO  173  K=l,NSIAR 
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z=o.o 

D3  150  J=l,NPl 
IND=K- Jf I 

IF  ( I.ND.LE.OI  GO  TO  155 
Z = Z«-B(  J)  ♦AIMAGI  YNYQAI  INO)  » 

150  CONTINUE 
155  00  160  J=l,NP 

IND=K- J 

IF ( INO.LE. 01  GO  TO  165 
Z = Z-A(  J»  ♦REAL  ( YNYQAI  INCH 
160  CONTINUE 
165  X=AIMAG( YNYQAlKn 

yNyOA(K)=CKPLX(Z,  <» 

170  CONTINUE 

00  180  K*l,NSTAR 
X=REAL  lYNYOAIKH 
130  YNYOAIK) =CMPLX( X,0.0t 

PRINT  OUT  THE  FILTERED  TIME  SERIES 

NS2=2^NSTAR 
WRITEINWf 12) 

)*R  I TEINW,  13)  (ARRAY  I I ) , I = I , NS2 , 2 ) 

12  F0RMATI30X,*  FILTERED  TIME  SERIES',/) 

13  FORMAT  ( 8( IXfE 13. 7) ) 

RETURN 

END 
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CnvPILER  OPTIONS  - NAME=  MA I N , OP T= 02 , L I NECNT =5 7 t S I ZE=OOOOK , 

SOURCE ,EBC01C ,NUL 1ST ,OECK,LOAD,MAP,NUEDIT,  ID,N0XREF 
SUBROUTINE  LNYQ 
COMMON  YNYQA(2043I 

COMMON  NW,NR , XF (2  0431  , YF 12  048  I , NS  I AR , NPOW , ND , JF , C T , ST , XOt  Y0» 

1 <S,YSfT,V,XNYO 
2fNT,XZ,YZ,  Jl) 

C JMMON  INTER PrlFFFT, IFFWAL , IFAUTO, IFGOF, IFPUNF,  IFPUNW, I FPUNA, 
I IFPUNO, IFPUNI , IFPUNR 
COMPLEX  YNYQA 

niMCNSICN  XOUT (4» , Y0UTI4)  x 
C THIS  ROUTINE  INTERPOLATES  AT  THE  NYQUIST  INTERVALS 
C USING  A STRAIGHT  LINE  AS  THE  INTERPOLATING  POLYNOMIAL 
C IF  XF I I-ll<X<XF(  I I THEN  THE  INTERPOLATED  VALUE 
C IS  Y V^HERE  Y = Cl*X*C2  WITH 
C Cl  = l YPl  n-YFI  I-m/KFlI  |-XF(  I-l)  I 
C C2=(Y=  (1-1  I^XFI  I )-YF(  I |*XFI  I-m/(XF  m-XF(  I-l)  ) 

IFLAG=-l 
L=  I 
1 = I 

20  IFIXFI I ) .GE.OIGO  TO  30 

1=  H-l 

IF ( I-JF) 20f 15C, 150 
30  JST=I 

YNYOAI L) =CMPLX(0.0,YF( I I ) 

4 0 <=  XF  I JST  )«-L»<NfQ 

90  I P I X-XF ( I ) ) IOC  , I 30, 140 

100  IFLAG=IFLAGU 

I 10  DENaM=  XF  I n-XF  I I-l  ) 

Cl=(YF(I )-YF( l-l) I/DENOM 

C2  = ( YF ( 1-  1) *XF I I ) -YF ( I l*XF(  I-l  ) ) /DENOM 

115  Y=Cl*X*C2 

116  l=L«-l 

IF (L-NSTAR  >120,120,  150 
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120  YNY0AlL)=CKPLX(0.0, Y» 

GH  TO  43 
130  Y=  YF ( I ) 

G3  TO  116 
143  IFLAG=-l 
1 = U I 

IF  I 1-JF I 90,90t 145 
145  l = JF 

GO  TO  no 
150  Wi^ITE(NW,3» 

WP I TE( NW  . 1 I 
30  163  K=l,NSTAR,4 
00  155  J=l,4 
JP  = 

X01IT(  J)  = XF(  JSTM-(  JP-ll*XNYQ 
155  YOUT(J)=AIMAG(YNYOA( JP)» 

WRITEINW,2M (XOUTlLI ,Y0UTCLI »,L=l,4» 

163  CONTINUE 

2 F0RMAT(3(2X,E13.7I» 

I FORMAT (4 ( I I X, • TIME  • ,6X  , • INTERP.  Y'l) 

3 FORMAT ( IHl ,20X, 'L ISTING  OF  THE  LINEARLY  INTERPOLATED  TIME  SERIES'! 
RETURN 

END 
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ns/360  FORTRAN  H 


COMPJIER  OPTIONS  - NAME=  MA I N , OP T = 3 2 t LI  NEC NT  = 5 7 » S I 7 E = 00D3K . 

SOURCE ,EBCDIC,NOL 1ST ,OECK,LOAO,MAP,NOEDIT, IO,NO<REF 
SUBROUTINE  NYQST 

C NVOST  TAKES  NON-EOUI  SPACED  ARRAY  YF.XF  AND  PRODUCES  AN  EQUISPACtO  COMPLEX 
C ARRAY  YNYQA.  INITIALLY  REIYNYQAI  =0  AND  I M I YN YOA ) = I NT E RPOL A TED  VALUES 
C EQUISPACEI)  ARRAY. 

C YFtn,  1 = 1, JF  GIVEN  AS  V/DIV 
C XFIII  , 1=1, JF  GIVEN  AS  NS/DIV 

C XNYQ  IS  the  NYOUIST  TIME  INTERVAL  IN  NS/DiV 
C NSTAR  = 2**NPQW  IS  GIVEN 
C 

C INITIALIZE  VALUES 
C 

DIMENSION  <DUT  K)  , YCUTI  A > 

COMMON  YNYCAI2DA8) 

CDMMON  Nri,NR,XF I20A3I ,YF(20A3» , NS T AR ,NPOW , NO , JF , C 1 , S T , <0, YO , 

I XS,YS, T, V,XNYQ 
2,NT,XZ,YZ, JO 

COMMON  INTERP,IFFFT,IFFWAL,IFAUTU, IFGOF , IFPUNF, IFPJNh , I F PUN A, 

1 1 FPUND, I FPUNI , I EPUNR 
COMPLEX  YNYQA 
DIMENSION  AI A,AI ,fiIAI 
L TSAV  = 0 
I 3 SAV=D 
I SAV=0 
I = I 
I ST*l 
C 

C SEARCH  FOR  FIRST  XF  > 0 

C 

20  IFIXFI  ISn.GE.OIGO  TO  30 

I ST=ISTf 1 

IF  I 1ST  - JFI  20,160,150 
30  JSTRT=IST 
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YNYQACLI  =CMPLX(O.Of  YF(  JSTRD  ) 

1= JSTRT 

70  XVAL  = XF(  JSTRT»  L*XNY0 
75  IF (XF( II .Gr.XVALIGO  TO  33 

1 = H-l 

IF  II  - JF)  75,115,75 
83  LT=Ul 

LB=I-2 

IFILB  .LE.  JSfRTIGO  TO  I 10 
1 = 1-1 
85  L=LM 

C 

C CHECK  IF  BOUNDS  ARE  SAME  AS  LAS!  TIME 

C 

IF (LT.E3.L rSAV  .AND.  LB.EQ.LBSAV  .AND,  I.EO.ISAVIGO  TO  100 
C 

C NEW  BOUNDS 

C 

LTSAV=LT 
lriSAV=L3 
I SAV=I 

IFIL  - NSTAR)  95,95,123 
C 

C FIT  THE  FOUR  POINTS  TO  A CUBIC 

C 

95  CALL  LSOPOL(A,B,2,LB,LTI 

DENOM=A( I , 1) ♦Al 2 ,2 I-A ( I, 21 *^2 
Tl=AI2,21*8m-A(  l,2M'B(2) 

T2  = A(  I ,21*81  I )-A( l,l )*BI2) 

B(  l)  = Tl/l)£NOM 
BI 2)=-T2/DENOM 

100  YNYQA(Ll=CMPLX|0.0,P0LYlb,2,XVAL) I 

IF (L-NSTARI 70,120,120 
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C NEAR  BEGINNING  OF  TRACE 

C 

llO  LR=JSTRT 
LT=L3f 3 
I=LB 

GD  TO  85 
C 

C NEAR  ENO  OF  TRACE 

C 

115  L9=JF-3 

L T = JF 
1 = I-l 
GO  TO  85 

150  WRITE(NM,II 

1 FORMATILX,'  ALL  TIMES  ARE  NEGATIVE') 

STOP 

120  WRITE(NW,3I 
WR  I TE  (NV,  ,4) 

00  200  K=l,NSTAR,4 
DO  190  J=l,^ 

INO=K^ J- I 

XOUTI  J ) = XF  ( JSTRTI  I I NO- I ) *XNyQ 
190  YOUn  J )=AI  MAGI  YNYQAI  IND)  I 

WRITE(NW,2M  ( (XQUriL)  ,Y0UTIL)  ) ,L=l,A)  ) 

200  CONTINUE 

2 F0RMAT(3I2X,F13.7) I 

3 FORMATIlHl , 'LISTING  OF  LEAST  SQUARES  INTERPOLATED  TIME  SERIES') 

A F ORMAT I A( I IX, ' TI ME ' ,6X, ' I NTERP.  Y ' ) ) 

RETURN 

END 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  M A I N , OP T = 02 i L I NECNT =5 7 , S I Z E = OOOOK , 

SOURCE.EBCOICfNOLIST ,OECK, LOAD, MAP, NOEOIT, ID,NOXREF 
SUBROUTINE  AUTCOR ( K ,NSAVE ,NPSAV) 

COMMON  YNV0A(2043» 

COMMON  NW,NR,XFI2  0't8I  ,yF(2048l  ,NS  T AR  , NPOW,  NO,  JF , C T , ST  , XO,  YO , 

I XS,YS, T, V, XN YO 
2, NT  ,XZ  ,YZ,  JD 

C 3MMON  INTERP, IFFFT , IFFWAL, I FAUTO,  IFGOF, IFPUNF, IFPJNW, IFPUNA, 
IIFPUND,IFPUNI,IFPUNR 
COMPLEX  YNYQA 

DIMENSION  FREQ (20481 , POWER (2348  I , ZF ( 2348  I , OUT ( 8 I 

DIMENSION  ARPAY(4096I 

EQUIVALENCE  ( ARRAY ( I I , YNYQA (II  I 

EQUIVALENCE  ( F R E 0 ( II  , ZF ( I II  , ( POWER! I I ,XF ( 1 1 ) 

WRITE(NW,l» 

I FORMATCl  LISTING  OF  THE  AUTOCORRELATION  FUNCTION* I 

00  PI=3. 1415927 

COMPUTE  AND  PRINT  OUT  THE  VALUES  OF  NSAVE,M,AND  IH 

K M l = K-  I 
IH=l 

M=NSAVE-NSAVE*KM1/  10 
NS  TAR  = M 

WRITE(NW,lOOl)  NSAVE,IH,M 
COMPUTE  DTAU^IH^XNYQ 
DT AU=IH*XNYO 

1001  FDRvaTI*  THERE  ARE*, 13,'  VALUES  IN  THE  DIGITAL  TIME  SERIES*,/, 

I*  THE  LAG  INTERVAL  IS  DTAU=IH*DT  WHERE  I H= * , I 3, / , * THERE  ARE*,  13,* 
2VALUES  IN  THE  DIGITAL  AUTOCORRELATION  FUNCTION*,/! 

C 

C COMPUTE  THE  AU TOCORREL A T I ON  FUNCTION  AND  STORE  IT  IN  THE  ARRAY  XF 
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c 

03  30  IR=l,M 
SUM=0.0 

I NMHR  = NSAVE  - IH*(  [R-U 
DO  20  IQ=l,INMHR 
IMESS=  IQ«-IH*(IR-l  ) 

20  SJM=SJMf YF 1 1 0) ♦YE ( I MESSI 

SUM=SUM/ INKHR 
XF( 1R»=SUM 
30  CONTINUE 

C 

C PRINT  OUT  AND  PLOT  THE  AUTOCORRELATION  FUNCTION 
C 

5 

VO 

VO  35 


36 

69 

A 

37 


C 

C COMPUTE  THE  RAW  SPECTRAL  DENSITY  ESTIMATES  AND  STORE  THEM 
C IN  THE  ARRAY  ZF. 

C 


WRITE(NW,5» 

FORMAT  I I X, AI lOX, • TI ME • ,7X, • AUTOCOR* ) I 
DO  35  1 = 1, M 

YNYQAIII =CMPLXl XF 111,0.01 

XNYSV=XNYQ 

XNYQ=DTAU 

DO  37  1=1, M, A 

DO  36  J=l,A 

OUTl 2* J-l I =XNYQ«l I -I fj-l ) 

□UT(2*JI=XF( I+J-l I 
WRITE(NW,AI  ICUTI LI ,L=l,8) 

FORMATIAI 7X,E13.7) I 
FORMATI IX,8I IX,E13.7I I 
CONTINUE 
CALL  PLOTSIAI 
XNYQ=XNYSV 


00  50  IR=l,M 
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is 

IPO=IR-l 

D 

SUM=XF(  UfXF(M)*(-l.»**IPO 

H 

X 

MMl=M-l 

DD  40  1Q=2,M«1 

> 

40 

sjM=siJM*2.o*xF  no)*cosiPi*(  [o-n*(iR-n/Mj 

l~  ( IR»=3rAU*S'JM 

50 

r 

CONTINUE 

L 

c 

COMPUTE  THE  REFINED  SPECTRAL  DENSITIES  AND  STORE  THEM  IN  THE 

c 

c 

ARRAY  Xr 

XF  ( 1 ) = 0. 54*2F  I lH-0. 23*ZFI2I 
XF  ( H » = 0.  23»ZF(  M-l  IfO.  54*ZF  (M» 
03  63  1=2, KMl 

60 

XF  m = 0.  2 3*ZF(  I-l  » f0.54*ZF(  I ) f 0.  23*ZF  ( I f 1 1 

03  70  1=1, M 

O 

o 

73 

r 

ZF  ( I I = (I-ll/(2*M*DTAU» 

v» 

C 

r 

PRINT  OUT  AND  PLOT  THE  REFINED  SPECTRAL  DENSITIES 

WRITE(NW,7» 

7 

f^ORMAT  ( • IREFINEO  SPECTRAL  DENSITT  ESTIMATES  FOR  THE  AUTOCORRELA T 10 
IN  FUNCTION'! 

RiRl  TE(NW,3I 

3 

FORMAT ( IH/ ,41 5X, ' FREQUENCY • ,9X , • POWER' 1 , /! 
DO  80  I=I,M,4 

80 

KR I TE(NW,4I  IFREQI I t J- 1 ! , POWE R ( I fj-l ) , J=l ,4! 

NSTAR=2*M 

<N  YSV=XNYO 

XNY0=3TAU 

00  93  1=1, M 

90 

YMYQAI I ) =CMPLXI POWER! I ! ,0.0! 
CALL  PL0TS(2! 

XN YQ=XNYSV 
NS  TAR=NSAVE 
RETURN 
END 

' ( JAN  73  I 


OS/363  FORTRAN  H 
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COMPILER  OPTIONS  - NAME=  MA I N , 0PT=02 t L I NECNT =5 7 t S I Z E=0000K i 

SOURCE .EBCOIC.NOL 1ST ,OECK,LOAD,MAP,NOEOITf ID.NOXREF 
SUBROUTINE  FFT 
C 

C THIS  ROUTINE  CALCULATES  THE  FOURIER  TRANSFORM  OF  A NYQUIST  SAMPLED  ARRAY 
C THAT  HAS  ITS  ELEMENTS  IN  BIT  REVERSED  ORDER.  THE  TECHNIQUE  IS  GENERALLY 
C KNOWN  AS  THE  FAST  FOURIER  TRANSFORM. 

C 

COMMON  YNYQA(20A3) 

COMMON  NW,NR,XFI2048I ,YF(2048I , NS T AR , NPOW , ND , JF, C T , ST , XOt YO , 

I XS.YS, TtV.XNYQ 
2,NT,XZ,YZ, JD 

COMMON  INTERP.IFFFT, I FF W AL , I F AU TO t I FGOF , IFPUNF, I F PUNW , I FPUNA , 
IIFPUNO, IFPUNI, IFPUNR 
COMPLE  X YNYQA, W, WE , F AC TOR , S WAP , 2 I P , S 
TO  I =6.  2831854 
DO  40  I=l,NPOW 
TI=2*»I 
ARG=TPI/TI 

S = CMPLX(COSI  ARGI  , SIM  ARGI  I 
IMI=I-l 
I l=2**IMl 
<STUP=I  If  I 

I 

W=CMPLX( 1 .0,0.0» 

20  SWAP=yNYOA (K) 

KO  I l=KH  I 

FACTOR  = YN/QA  I KPI  U«W- 
YN YQAlKI=SWAPfFACTOR 
YNYOAIKPI I )=SWAP-FACTOR 
K = Kf  I 
W = W*S 

IF IK-KSTQP)20, 30,20 
30  K=KPIlfl 
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102 


KST0P=Kf I I 
W=CMPLX( 1.0,0.01 
IF ( KSTOP.LE.NSTARJ  GO  TO  20 
40  CONTINUE 

DELF=l  ./(NSTAR*XNYQ) 

WR I TE( NW,  I I 

XS0=XNYQ**2 

DO  50  I=1,NSTAR,ND 

Z IP  = YNYQA( I I 

RZ IP=RFAL( ZIPI 

AI ZIP=AIMAGl Z IPI 

FREQ=( I-l»*OELF 

P/#R=(RZIP**2«-AIZ  IP**2»*XS0 

PHS=ATAN(AIZIP/RZIP) 

WR I TE(NW,2»FREQ,RZIP, AIZIP,PWR,PHS 
YNYQAl I |=CKPLX(PWR,PHS) 

50  CONTINUE 

1 FOPMATIlHl,'  FREQUENCY  REAL  F(W)  IMAG  FIWI 

lER  PHASE*) 

2 F0RMAT(5(IX,E13.7)) 

R5  TURN 

END 


POW 
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OS/360  FORTRAN  H 


COMPlLtR  OPTIONS  - NAME=  MA I N, OP T = 02 t L I NECNT  = 57 , S I 7 E = a330K , 

SOURCE ♦ EBCDIC »NOL I S T , DECK ,L0 AD » MAP, NOEO I T , ID,NOXRFF 
SUBROUTINE  FWAL 

THIS  ROUTINE  CALCULATES  WALSH  TRANSFORMS  OF  A NYQUIST 

SAMPLED  ARRAY.  THE  TECHNIQUE  IS  GENERALLY  KNOWN  AS  THE  FAST  WALSH  TRANSFORM. 
COMMON  YNYQA(2348> 

COMMON  NW,NR,<FI2  34B)  , T’F  (2  348)  , NS  T AR  , NPOW , ND , JF  , CT  , S T , X3,  YO, 

I XS,YS,T,V,XNYO 
2,NI,XZ,YZ,JD 

COMMON  INTERP,  IFFFT, IFFWAL , IFAUTO, IFGOF, IFPUNF, IFPUNW, IFPUNA, 
IIFPUNO.IFPUNI , IFPUNR 
COMPLEX  YNYQA 
DO  10  I=l,NSTAR 
XF I I )=REAL( YNYQAI  I ) ) 

10  YF(I)=0.0 

DO  40  1=1, NPOW 
IM  l=I-l 
I l = 2»*IMl 
KST0P=1  U I 
L=  I 
K=  I 
J = 0 
W=l.3 

23  Sl=<F(K) 

KP  I l = KH  I 
S2  = XF( KP  I 1 ) 

SW2=S2*W 
INDEX= J/2 
SIG=(-ll**INDEX 
YF(L)=(SUSW2)+SIG 
YF  (LH  ) = l-SH-SW2)  *SIG 
K = K>  I 
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J=  J*^l 
L = Lf2 
W=-W 

IF(K-KST0P)20,30,20 
30  K = KPIUl 

KSTOP=K«-I  I 

J = 0 

K=l 

IF(KSTOP.LE.NSTAR»  GO  TO  20 
00  35  IZX=UNSTAR 
XF(IZX»  = YFnZX» 

35  VFnZXI=0.0 

40  CONTINUE 

DELF= I , / INSTARtXNYQJ 
WR ITEINUt  1) 

WR ITE(NWt2) 

NSMl=NSTAR-l 
G=XFI l»**2 
L = 0 

YNYQA(Ltl)=CMPL<lG,0.0» 

WR(TE(NW.3l  LfXFm.G 
00  50  I=2,NSMl,2 
L=L*l 
AS  = XF( I I 
AC=XFI 1*1) 

G=AS**2«-AC**2 

YNYQAI  LH  » =CMPLX(G»0.0» 

50  WRITE(NWf4)  L»ASfAC»G 

L = L*l 

G=XFINSTAR)**2 
YNYOAILUI  =CMPLX(Gf0.0) 

WRITE(NW»5I  L f XF(NSTAR) «G 

1 FORMATIlHl,'  OUTPUT  OF  THE  WALSH  TRANSFORM’ » 

2 FORMATI’  I S9X, 'ASin  »»9X,  »ACm  *,10X, 'Gm*  » 


> 

v 

V 


a 

H 


X 


> 


rf\  -t  \r\ 


F0RMAT(K,I5,14X,2(  I <♦  El  3. 7 II 

FORMAT(lX»I5»3llX,El3.7l I 

F0RMATllX,I5»lXtEl3.7tl5XtEl3.7l 

RETURN 

END 


c 

Lr. 
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OS/360  FORTRAN  H 


COMPILER 

OPTIONS  - NAME=  MA I N t OPT^OZ » LI N£CNT= 5 7, S I ZE=0000K , 

S0URCE,EBC0IC,N0L 1ST, DECK, LOAD, MAP, NOEOIT, I0,N0XREF 
SJBRCUTINE  INVMATI A,6,N,PIV0T, INV) 

01  ME  NS  ION  AK,4I  ,B|4I  , INDEX(^,2)  , IPIVOT  (4) 

EQUIVALENCE  I JROW,IROW), 1 JCOLUM, I COLUM » , 1 AM AX,  SWAP, T I 

C 

c 

c 

c 

c 

c 

c 

MATRIX  INVERSION  WITH  ACCOMPANTING  SOLUTION  OF  LINEAR  EQUATIONS 

INITIALIZATION 

DO  10  1=1, N 

IPI VOT ( I) =0 

10 

♦-> 

CONTINUE 
00  100  1=1, N 

90 

o o o 

SEARCH  FOR  PIVOT  ELEMENT 
AMAX=0.0 

20 

DO  40  J=l ,N 

IFIIPIVOTUI  .EQ.U  GO  TO  40 
00  30  K=l,N 

IFI  IPIVOTIKI-I  120,  30,  130 

IF( ABSI AMAXI .GT.ABSI A( J,K» ) I GO  TO  30 

IROW=J 

ICOLUM=K 

AMAX=A( J,K) 

o o 

CONTINUE 

CONTINUE 

IPI  VOT  ( I COLUM  ) = IPI  VOT  ( ICOLUMM-l 

c 

c 

c 

INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

APPENDIX  A 


o o o o o o 


[F ( IROW.EQ. ICOLUM)  GO  TO  60 
00  50  J=l,N 

SWAP=A(IROW, J) 

AdROW,  J)  = A(  ICOLUM,  J| 

AdCOLUM,  J»  = SWAP 
50  CONTINUE 

SWAP=B( IROW ) 

B( IR0WI=B( ICOLUMJ 
B ( ICOLUM»=SWAP 
PI VOT=A  dCOLUM, ICOLUH» 

INDEX! I ,2»=IC0LUM 
60  INOEXCI ,l»=IROW 

IF(PI VOT.EO.O.OI  GO  TO  130 

DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

A! tCOLUM, ICOLUMI=l.O 
DO  70  J=l,N 

AdCOLUM,  J»=A(  ICOLUM,  JI/PIVOT 
CONTINUE 

3 ( IC0LUM»=B( ICOLUMI /PIVOT 
REDUCE  NON-PIVOT  ROWS 
DO  90  J=l,N 

IF ( J, EQ. ICOLUMI  GO  TO  90 
T=A( J, ICOLUM) 

A! J,1C0LUMI=0.0 
DO  30  K=I,N 

A(J,K»*A(J,K»-A( ICOLUM,KI*T 
80  CONTINUE 

B(J»=B( J)-B( ICOLUM)*T 
90  CONTINUE 

100  CONTINUE 
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non 


IV1T6RCHANGE  CCLUMNS 


o 

00 


I 10 
IdO 
I 30 


DO  120 


RE  ru«M 
FS!0 


I =l  »N 
J=NH-1 

IF( INOE <( J , n .EQ. INDEKl J, 2 ) ) GO  fO  120 
JttQW=INDEX (J, n 
JC0LUM=IN0EX(J,2I 
•)(1  HO  K=l,N 

SWAP=A( K, JROWJ 
A(K,  JR0WJ=A(K,JC0LUM1 
AIK,  JCOLUM)=SWAP 
CONT INUE 
CONTI NLE 
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OS/360  FORTRAN  H 


OMPILER  OPTIONS  - NAME=  M A I N , OP T*02 i L I NEC NT  = 5 7 , S I ZE=OOOOK , 

SOURCE .EBCDIC »NOL 1ST .DECK , LOAD. MAP, NOED I T , ID.NOXREF 
SUBROUTINE  L SCPOLI A . B . N , L B , L T » 

C 

C THIS  ROUTINE  CALCULATES  THE  MATRIX  A AND  VECTOR  B WHICH  ARE  USED  TO  SOLVE  FOR 
C THE  VECTOR  X IN  THE  EQUATION  A,X=8.  THE  VECTOR  X HAS  ELEMENTS  WHICH  CORRESPON 
C TO  THC  POLYNOMIAL  COEFFICIENTS  IN  SUMCOE FF I II ♦ Z* ♦ I . T HE  COEFFICIENTS  ARE 
C DETERMINED  IN  A LEAST  SQUARES  SOLUTION  OF  THE  POLYNOMIAL  TO  SOME  DATA, 

C 

common  YNf0A(20A8» 

COMMON  NW.NR  ,XF(  20<»9)  ,YF  ( 20481  . NST  AR  , NPQW  , NO  , JF  , C T , ST  , X 0 , YD, 

I X5.YS, r, V.XNYQ 
2,NT,XZ.YZ,  JD 
COMPLEX  YNYQA 

DIMENSION  A(4,4I ,B(4) ,C( 10) 

LQ=2*N-l 
NPTS=L  T-LBt I 
DO  13  1=1, N 
31  I »=0.0 
DO  10  J=I ,N 
ID  A(  I,  Jl=3.0 
DO  5 L = l ,L0 
5 C(U=D.D 

C(  l)=NPTS 
DO  15  l=Lfl,LT 
00  13  K=l, N 
KM l=K-  I 

IF  I XF(  I ) .EO.O.O.AND.KMl.EQ.O)  GO  TO  113 
BIKI=8(KM-YF(I  )♦(  XF(I  )*»KM1) 

GO  TO  13  > 

I 13  BI  K I =B  (K  ) f YF  ( II  ^ 

13  CONTINUE  g 

DO  14  L=2.L0  g 


14  C (LI=C  (U  + XF  ( n«*LMl 

15  CONTINUE 

00  20  K=l,N 
K'l  l=K-  I 
PO  20  KP=K,N 
MK,KP)=CIKP«-KMII 
20  A( KP,K )= A(K,KP » 

PF  TURN 
END 


o 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  MA I N , OP T =02 # LI NECNT =5 7 , S I ZE=OOOOK , 

SOURCE  tEBCOICfNOL 1ST ,D£CK, LOAD, MAP, NOEO IT, IO,NOXREF 
SUBROUTINE  PLCTSINI 

C***  PLOT  ONLY  THE  LOG  OF  THE  POWER  FOR  THE  MOMENT 
COMMON  YNYQA(2043» 

COMMON  NW,NR ,XF (20481 , YF ( 2048 ) , NS T AR , NPOW , NO , JF , CT , ST , XO, YO, 

I XS,YS, T,Y,XNYQ 
2,NT,XZ ,YZ , JO 

COMMON  INTERP, IFFFT, IFFWAL , I FAUTO, IFGOF, IFPUNF, IFPUNW , I FPUNA, 
IIFPUNO, IFPUNI, IFPUNR 
COMPLEX  YNYOA 

DIMENSION  ARP (126  I , PHI  (1024) ,Y(  1024) 

DATA  BLANK, STAR, ZERO, PLUS, BAD/IH  , IH*  , IHO,  IH*-,  IHB/ 

NCHAN=0 

GO  TO  ( 100,200,300,400)  ,N 
100  DO  llO  I=l,NSrAR,JD 
NCHAN=NCHANf I 

PHI (NCHAN) =AIMAG( YNYQA(  I ) ) 
no  Y(  NCHAN)  =PHI  (NCHAN) 

G3  TO  500 
200  Nl=NSTAR/2 

on  210  l=l,Nl,ND 

NCHAN=NCHAM-  I 

PHI ( NCHAN )=REAL( YNYOA ( I) ) 

210  Y( NCHAN) =PHI (NCHAN) 

GD  TO  500 

300  DO  310  I=l,NSTAR,ND  . 

NCHAN=NCHANf  I 

PH  I ( NCHAN )=AI MAGI YNYOA ( I ) ) 

310  Y( NCHAN) =PHI (NCHAN) 

GO  TO  500 

400  no  410  I=l,NSTAR,JD 
NCHAN=NCHAN<- 1 
PH  I ( NCHAN) =REAL ( YNYQA ( I ) ) 
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413 
5 30 


10 

c *** 


23 

30 

C<=** 


40 


53 
6 0 
635 
61 

62 

70 


Y(  NCH&N')  =PHI  (NCHANI 
NPSTEP=l 

THIS  ROUTINE  DOGS  THE  PLOTTING  OF  THE  DATA  AND  THEORETICAL  FIT. 
WR  I TE  ( NVs,  n 

INITIALIZE  FOP  FINDING  MAXIMUM  AND  MINIMUM. 

4MAX  = PHI  I 1) 

AM IN=AMAX 

INITIALIZE  THE  PLOTTING  ARRAY  TO  BLANKS. 

01  13  1=1,126 

ARR(  n=BLANK 
CONTINUE 

FIND  MAXIMUM  AND  MINIMUM  OF  EITHER  THE  DATA  OR  THE  FIT. 

33  33  I=l,NCHAN,NPSTEP 

IF ( YI n .GT.AMAX)  AMAX=Y(II 
IF  ( Y(  I I .LT.AMINI  AMIN=Ym 
IF  (PHI  ( I I .GT  . AMAX)  AMA<=PHI(.I» 

IF (PHI  ( II  .LT.AMINI  AMIN=PHI(I) 

CONTINUE 

FIND  THE  BIN  SIZE. 

D3IN=( AMAX-AMINI/ 124.3 
FIND  WHERE  TO  PUT  THE  CURVES. 

DD  8C  I = I ,NCHAN,NPSTEP 

I I =IFI X( ( PHI ( I I -AMINI ZDBINI *2 
JJ=IFIX((Y(I)-AMINI ZOBIN) f2 
IF  ( I I.NE. JJ I GO  TO  50 
ARRd  I l=PLUS 
GO  TO  635 
ARR( JJI=ZERO 
ARP ( II l=STAR 
GO  TO  (6  I ,62,62,61  I ,N 
K=(  I-l )*JD 
GO  TO  73 
K=( I-l )*ND 

WRITE(Nls,2I  K,  (ARR(KI  ,K=l,l26) 
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ARRl  I I »=RL4NK 
ARR( JJI=BLANK 

80 

CONTINUE 

DELF=l,/(NSTAR»XNYOI 
G3  TO  (9lf82,83,34) ,N 

81 

WP  I TE(NW,4I 

4 

fOR'IATC  PLOT  OF  AMPLITUDE  VS  TIME  AFTER  INTERPOLAT  ION  • » 
WRITE(NW,5)  XNYQ 

5 

FORMAIC  TO  OBTAIN  TIME  VALUES  MULTPLY  XCOORD  BY*,E13.7I 
RE  TURN 

8^ 

WR I TE(NW,6> 

6 

FORMAT!'  PLOT  OF  THF  FILTERED  TIME  SERIES'! 

W«ITE(NW,5I  XNYO 

RETURN 

92 

WR I TEINK, 7» 

M 

7 

FORMAT!'  PLOT  OF  THE  POWER  VS.  FREQUENCY'! 

H 

OJ 

83 

WRITE!NW,3!  DELF 

RETURN 

WR I TE!NK,8  ! 

8 

FORMAT!'  PLOT  OF  THE  PHASE  VS.  FREQUENCY'! 

WRITE(NW,3!  DELF 

RETURN 

3 

FORMAT!'  TO  OBTAIN  FREQUENCY  VALUES  MULTIPLY  X VALUE  BY',E13.7! 

1 

FORMAT  I IHl! 

2 

FORMAT  ! IX, I 4, 2X,  126AI ! 
END 

i§ 

n 

z 

3 

H 

X 

\ 

V, 

\ 

1 • 

> 
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OS/363  FORTRAN  H 


compiler 

OPTIONS  - NAM6=  M A I N , OP T = 0 2 » L I NECNT  = 5 7 , S 1 ZE=0000< , 

SOURCE , EBCDIC ,N0L 1ST , DECK, LOAD, MAP,N0EDI T, ID,N0XREF 
SUBROUTINE  CSTOUT 

T 

o o o 

CHECK  TIME  ORDERING  OF  THE  POINTS  AND  CAST  OUT  THOSE  POINTS  NOT  IN  THE 
PROPER  TIME  ORDER 

IT  4LS0  AVERAGES  THOSE  POINTS  WHICH  HAVE  THE  SAME  < VALUES 
COMMON  YNYQA(20^8I 

COMMON  NW,NR,XF (2048) ,YF 12048 ) ,NS T AR , NPOW , ND , JF , CT , S T , XO, YO, 

1 XS,YS,T,V,XNYO 
2,NT,XZ,YZ, JD 

COMMON  INTERP,IFFFT, I FF WAL , I F AUT 0, I FGOF , IFPUNF, I FPUNW , I FPUNA , 

IIi^PUNO,  IFPUNI  , IFPUNR 
COMPLEX  YNYOA 
LSuy=i 
J=l 

J0AD=3 

K=2 

YSUM  = YF I 1 ) 
XSAV=XF( I ) 

80 

90 

IF  1 XFI K 1 -XSAV) 90, 100, llO 
K = K + l 

JBA()=JBAD  ♦ 1 
IF(K-JF)83,33,113 

100 

LSUM=LSUMfl 
YSUM=YSUMf YF IK ) 

K = K U 

IFIK-JF)80,80,30 

1 10 

YF  I J ) = YSUM/L  SUM 
XFI J)=XSAV 
IFIK-JFI 120, 133, 33 

120 

XSAV=XF( K) 
YSUM=YF ( K ) 
l.SU'^=I 
K = K ♦ 1 
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J-  J»- 1 

DU  TO  80 
130  J = JH 

YF  ( J)  = YF  (K> 

XF  ( J|  = XF (K) 

30  WPITE(NW,l)  JBAD,JF»J 

I FORMAT (/ , IXt I 5 , • POINTS  OUT  OF  A TOTAL  OF  • t 15, 

1 • on  NOT  SATISFY  THE  TIME  ORDER  CRITERION*, 

2 /,I5,'  POINTS  WILL  BE  USED') 

JF=J 

RETURN 

END 


IT. 


X 

> 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  MA I N ,0PT»02, L I NECNT=5 7 , S I ZE=OOOOK , 

SOURCE, EBCO I C,NOL 1ST ,DEC K, LOAD, MAP, NOED I T , ID,NOXREF 
SUBROUTINE  SAVE  INI 
COMMON  YNVQA (20481 

COMMON  MW,NR,XF(2048» , YF I 2048 » , NST AR , NPOW ,ND , JF , CT , ST , XO, VD, 

I XS,YS,T,V,XNYO 
2,NIT,XZ,YZ,  JO 

COMMON  I NTERP, I FFF T , I FF WAL , I FAUTO , I FGOF , IFPUNF , I FPUNW , I FPUNA, 
IIFPUNO, IFPUNI , IFPUNR 
COMPLEX  YNYQA 
IF(N)  10,10,20 
10  DO  15  I=1,NSTAR 

15  YFm^AIMAG  (YNYQAmi 

RETURN 

20  00  25  I=1,NSTAR 

25  YF(I»»REAL  (YNYOAIIM 

RETURN 
END 
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OS/363  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  M A I N , OP T*0 2 , L I NECNT  = 57 , S I ZE= OOOOK , 

SOURCE , EBCDIC »NOL 1ST , DECK, LOAD, MAP, NOEO IT, ID,NOXREF 
SUBROUTINE  UNSAVEIN) 

COMMON  YNfQAI2043l 

COMMON  NW, NR ,XF (20481  , YF I 2048 » , NS T AR , NPOW , NO , JF , C T , ST , XO , YO, 

I XS, YS,T,V,XNYQ 
2,NT,KZ, YZ, JD 

COMMON  I NTERP, IFFFT, I FF WAL , I FAUTO , I FGOF , I FPUNF, IFPUNW, IFPUNA, 
IIFPUND,IFPUNI, IFPUNR 
COMPLEX  YNYQA 
IF  (N)  10,10,20 
10  DO  15  I=l,NSTAR 

15  YNYQAm=CMPLXlO.O,YFim 

RETURN 

20  DO  25  I*l,NSTAR 

25  YMYOAI  n=CKPLXlYF(  11,0.01 

RE  TURN 
END 
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( JAN  73  I OS/360  FORTRAN  H 

COMPILER  OPTIONS  - NAME=  MA I N » OP T=02 , L I NECNT =5 7 , S I ZE= 0300K , 

SOURCE  tEBCOICtNOL 1ST ,DECK,LOADf  MAP, NOEDIT*  IO,NOXREF 
SUBROUTINE  ROT 

C ROTATE  ZERO  POINT  OF  THE  TRACE 

COMMON  YNYQAI20A3) 

COMMON  NW,NR,XF (20A8>  , YF( 2048  I , NS T AR , NPOW , NO , JF , C T , ST , XO , YO, 

I XS,YS,T,V,X\YO 
2,NT,XZ ,YZ , JO 

COMMON  INTERP, IFFFT , IFFWAL,IF AUTO, IFGOF, IFPUNF, IFPUNW, IFPUNA, 
IIFPUNO, I FPUN I , IFPUNR 
COMPLEX  YNYOA 
XZS=XZ-XO 
YZS=YZ-YO 

XZ=  CT#XZS  *■  ST*YZS 
YZ=  -ST*XZS  ♦ CT*YZS 
C ROTATE  AND  SCALE  INPUT  ARRAY 

DO  10  1 = 1, JF 
SW<=XF ( I ) - XO 
SWY=  YF ( n - YO 

YFm=  ( -ST*SWX*-CT*SWY  - YZI*V 
10  XF  ( n = (CT*SWX+ST*SWY-XZ)  *T 
RETURN 
END 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME-  MA I N, OPT *02 , L I NECNT-5 7, S I ZE-OOOOK , 

SOURCE fEBCOIC  »NOL 1ST ,OECK, LOAD, MAP,NOEOIT, IO,NOXREF 
SUBROUTINE  SCARTP 
COMMON  TNTQA(2048I 

COMMON  NW, NR ,XF (20481 , YF ( 2048 » , NST AR , NPOW, NO, JF ,CT , ST , XO, YO, 

1  XS,  YS,  T,  V,  XNYQ 
2,NT,XZ,YZ, JO 

COMMON  INTERP, IFFFT, I FFW AL , I FAUTO , IFGOF , IFPUNF, IFPUNW , I FPUNA, 
lIFPUSOtlFPUNI, IFPUNR 
COMPLEX  YNYOA 
DIMENSION  XT(4»,YT(4» 

PEAD(NR,2)XZ,YZ 
READINR,U  Xl,X2,Yl,Y2 

1 F0RMATI4I5I 

REAOINRfZ)  (XTII  I ,YT( I), 1 = 1,4) 

2 F0RMAT(2( IX,E11.4) ) 

RE A0(NR,2)T,V 

3 FORMAT! IX, • THETA*  * ,E13.7/ IX , 'COS  I THE T A I = * ,E13.7/ 

X IX, 'SINITHETAI-SEIJ.T/IX, ’COORDINATE  ORIGIN  S = ',E13.7/ 

X I X, 'COORDINATE  ORIGIN  Y = ' , E 1 3 . 7, I X, • SC AL E FACTOR  X » »,E13.7/ 

X IX, 'SCALE  FACTOR  Y = ’,E13.7/I 

4 FORMAT  I I X, ' XAXI S PTS  X,Y  COORDINATES’,/, 

X 3( IX,E13. 7) , /, 3( IX,E13. 7) ) 

5 FORMAT  I IX,  *YAXI S PTS  X,Y  COORDINATES',/, 

X 3IIX,EI 3. 7) ,/,3( IX,E13.7) ) 

XO  = ABS  I X2*XT(  1 1 -Xl»XT  (2)  )/ABS(Yl-X2) 

YO-ABS  (Y2*YT(3)-Yl*YT(4»  )/ABS(Yl-Y2) 

XS-SORTI (XTI 2)-<TI  1) ) ♦♦2* I YT I 2 I - Y T I U l♦♦2)/ABS  (X2-XII 
YS»SQRT  I ( XTI4)  - XT  ( 3)  ) ♦♦2»- ( Y T ( 4 ) -Y  T I 3 I )**2»/ABS  I Y2-YI ) 

IF  I YT( 2) .EQ. YI I 1)  I GO  TO  20 

TANT=(  (YTI  2)-YT(l)  ) / (XTl  2)-XT(  II  )-IXT(4)-XT(3n/IYT(4)-YTI  3m/2.0 
CT=l.O/SQRT( l,Of TANT**2) 

ST=TANT*CT 
THETA=ATANI TANTI 
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GO  TO  30 
20  THfT&=0.0 

ST=0.0 
CT=1.0 
30  T»T/XS 

V= V^VS 

ITE(NWt3l  THETA,CT,ST,X0,Y0f XSt YS 
WRMEINW.AI  Xl,XT(  U,YT(2l,X2,XT(2»,Yn2» 
WRI1E(NW,5)  Yl  ,XT13I,YT(3) ,Y2tXT(4)f YT(A) 
RETURN 
END 
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QS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAM6=  MA I N , 0PT*02 t L I NECNT = 5 7 t S I ZE=0000K t 

S0URCE,E8CDIC  f NOLI ST, DECK, LOAD » MAP, NOEDIT, ID,NOXREF 
SJBROUTINE  SRTFUR 

ms  ROUTINE  PERFORMS  A BIT  REVERSAL  OF  THE  ORIGINAL  NYQUIST  ARRAY. 
IF  I IS  THE  INDEX  OF  AN  ELEMENT  OF  THE  ARRAY  THEN  I = SUM I A I N I ♦ 2**N I 
WHERE  THE  LIMITS  OF  THE  SUM  ARE  0 AND  M AND  AIN>  IS  0 OR  1. 

THE  ELEMENT  SPECIFIED  BY  I IS  EXCHANGED  WITH  THE  ELEMENT  SPECIFIED 
BY  J,  WHERE  J=SUM( A( M-N» ♦2*»N. 

EXAMPLE:  SUPPOSE  1=0110101,  THEN  J=lOlOllO. 

NS  TAR  = 2**NPOW 
COMMON  YNY0AI2048I 

COMMON  N W, NR ,XF (20481 ,YF ( 2048  I , NS T AR , NPOW , ND , JF , CT , S T , XO , YO , 

1 XS, YS, T, V, XNYQ 
2,NT,XZ,YZ, JO 

COMMON  INTERP,IFFFT, IFFWAL,IFAUTO, IFGOF, IFPUNF, I F PUNW , I FPUN A , 
IIFPUNO.IFPUNI, IFPUNR 
COMPLEX  YNYOA,SWAP 
DIMENSION  INT115I 
DO  10  I»l,NPOW 
10  INTI  U=2**INP0W-n 
DO  30  I=I,NSTAR 
J=  I-l 
I SUM=0 

DO  20  K=I,NP0W 
ISUM=ISUM  MODI  J,2I  ♦‘INTIKI 
20  J= J/2 

ISUM=ISUMf  I 
IF  { I SUM-I » 30,30,25 
25  SWAP= YNYOA I I I 

YNYOAII |=YNYOAI ISUM) 

YNYOAI ISUMJ =ShAP 
30  CONTINUE 

RE  TURN 
END 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  M A I N, OP T = 02 , L I NECNT  = 5 7 , S I ZE  = OOOOK , 

SOURCE  fEBCDIC  tNOLIST  ,()ECK,LOAD,  MAP,  NOEO  IT,  I 
FUNCTION  POLY(B,N,XVAL) 

EVALUATES  A POLYNOMIAL  FUNCTION 

DIMENSION  BU) 

POL  Y = 3.0 
XP=l .0 
DO  10  1=1, N 
POLY  = POLY  B( n »XP 

13  XP=XP*XVAL 
RETURN 
END 


H 

n; 

K; 


D,N0XR6F 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  MA I N t0PT»02 » LI NECNT*57 , S I ZE’OOODK * 

SOURCE, EBCOlCf NOLI  ST, DECK, LOAD, HAP, NOEDIT, ID,NOXREF 
SUBROUTINE  POLMLT ( M, A , N, B,L ,C I 
COMPLEX  CI40) ,B(20I , AI20I 
L=MfN-I 
DO  10  I*l,L 

10  C(  n=CMPLX(0.0,0.0» 

00  20  J=l,M 
DO  20  K=l,N 

1 = J*^K“I 

20  C(  II=C(I  l«^Al 

RETURN 
END 


APPENDIX  A 


APPENDIX  B.  SIGNAL  PROGRAM 


( JAN  73  ) OS/360  FORTRAN  H 

COMPILER  OPTIONS  - NAMB=  MA I N» OPT=02 t L I NECNT*5 7 » S I ZE=OOOOK , 

SOURCE  rEBCDIC, NOLIST, DECK, LOAD, MAP, NOEOIT, ID,NOXREF 
C 

C MAIN  PROGRAM  FOR  A GENERAL  NON-LINEAR  LEAST  SQUARES  FIT  OF  A 
C PARAMETRIZED  FUNCTION  TO  A SET  OF  DATA  Y 
C 
C 
C 

C NON  LINEAR  LEAST  SQUARES  FITTING  PROGRAM 

C 

C 

C 

C 

COMMON  IMFORMATI 

C A MATRIX  OF  EON.  8 AND  INVERSE  AFTER  INVERSION 

C C VECTOR  IN  EQN.8 

C DPR  CHANGES  IN  THE  RESTRICTED  PARAMETERS 

C FMT  FORMAT  FOR  THE  EXPERIMENTAL  DATA 

C PR  RESTRICTED  PARAMETERS 

C PSAVR  ARRAY  FOR  TEMPORARY  SAVING  OF  PR 

C Z . ARRAY  OF  DERIVATIVES 

C PHI  VALUE  OF  THEORETICAL  FUNCTION 

C RESID  RESIDUALS 

C X ARRAY  OF  INDEPENDENT  VARIABLES  SUCH  AS  TIME 

C Y THE  EXPERIMENTAL  DATA 

C W THE  STATISTICAL  WEIGHTS 

C IFCOR  DETERMINES  DISPOSITION  OF  CORRELATION  COEFFICIENTS 

C It^PLOT  DETERMINES  IF  THE  PLOTTING  PROGRAM  IS  CALLED 

C lYDEZ  INTERNAL  PARAMETER  COMMUNICATING  BETWEEN  MAIN  AND  TESTS 

C THE  INFORMATION  THAT  THE  PROGRAM  HAS  OR  HAS  NOT 

C CONVERGED 


o oooooooooooooonooooooooo 


N3  THE  NUMBER  OF  THE  BEGINNING  PARAMETER  IN  A GIVEN 

REGION  OF  PARAMETER  SPACE 
N3ADPT  THE  NUMBER  OF  BAD  DATA  POINTS 

N:HAN  the  NUMBER  OF  CHANNELS 

NE  THE  NUMBER  OF  THE  LAST  PARAMETER  IN  A GIVEN  REGION  OF 

PARAMETER  SPACE 

NP  NUMBER  OF  PARAMETERS  TO  BE  FIT  IN  A GIVEN  REGION  OF 

PARAMETER  SPACE 
N»ST£P  USED  FOR  PLOTTING 

VP  THE  NUMBER  OF  PARAMETERS  TO  BE  FIT  IN  A GINEN  REGION 

OF  PARAMETER  SPACE 

N^P  THE  NUMBER  OF  RESTRICTED  PARAMETERS 

AA  ALPHA  IN  EON.  12 

CHI  CHI-SQUARED 

CHIB  CHI-SOUARED  BEFORE  EACH  ITERATION 

CHIA  CHI-SQUARED  AFTER  EACH  ITERATION 

DEGF  THE  NUMBER  OF  DEGREES  OF  FREEDOM 

iP*Att**<f****fi^^  ******  ^iin^^nti*********************^**’^******************** 


♦♦♦  I/O  INFORMATION 

♦♦♦  NREAD  IS  THE  S/MBOLIC  DESIGNATION  FOR  THE  CARD  READER 

NORITE  IS  THE  SYMBOLIC  DESIGNATION  FOR  THE  LINE  PRINTER 
♦♦♦  NPUNCH  IS  THE  SYMBOLIC  DESIGNATION  FOR  THE  CARO  PUNCH 
***  i/Q  CHANNELS  ARE  DEFINED  IN  THE  NEXT  THREE  STATEMENTS 

CDMMON  AI50»50) ,C(50I ,0PRI50) ,FMTI5» ,PR ( 50 1 , PSAVR I 50  I f Z ( 50  I , P I 50» 
COMMON  PHI ( 1024) rRESI0(l024> »X( 1024) ,Y( 102^1 »H( 1024) 

COMMON  I FCORt IFPLOT,INOEZ(NB,NBADPT»NCHANtNE»NP»NPSTEP>NRP» 

1 NREADrNWRlTEiNPUNCH 

COMMON  AA,CHI »CHIB»CH1A»DEGF 
N3EAD=5 
NWRITE*6 
NPUNCH=7 
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GET  THE  MAXIMUM  # OF  ITERATIONS  AND  THE  INPUT  DATA  FORMAT 
READINREAOfll  JSTOP,FMT 
GET  THE  INPUT  DATA  FOR  THE  N-TH  FIT 
0 CALL  READIN 

INITIALIZE  ITERATION  COUNTER 
J0UIT=O 

WRITE(NMRITE»2» 

INCRE'IENT  ITERATION  COUNTER 

0 J3UIT=JQUIT+1 

TEST  n OF  ITERATIONS 

I*:!  JQUIT.LE.  JSTDPIGO  TO  30 
WRITE(NHRITE»3) 

GO  TO  ^0 

SET  UP  MATRIX  OF  DERIVATES  AND  SOLVE  FOR  DELTA  PARAMETERS 
0 CALL  GRINDIOI 

TEST  AND  MAKE  CHANGES  IN  THE  PARAMETERS 

CALL  TESTS 
C 

C IN0EZ=0  IF  THE  FIT  HAS  CONVERGED 
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c 

IF ( INDEZ ) 23,40,43 

WRITE  OUT  THE  PARAMETERS  AND  THE  THEORETICAL  FIT 
0 CALL  SCRIBE 

TEST  IF  A PLOT  IS  DESIRED 
IF(  IFPL0TI60, 70,60 
DO  THE  PLOTTING 
0 CALL  PLOTS 

DO  THE  ERROR  ANALYSIS 
0 CALL  ERRMAT 

GO  GET  THE  DATA  FOR  THE  NEXT  FIT 

IF (1)80,60,10 

1 FORMAT!! 5,5A4) 

2 FORMAT! IH-, 'CHANGES  IN  PARAMETERS') 

3 FORMATIIH-,'  THE  NUMBER  OF  ITERATIONS  EXCEEDED  JSTOP.  THE't/,' 

1 ' PROGRAM  WILL  AUTOMATICALLY  GO  TO  THE  PLOTTING  AND  ',/, 

2 • ERROR  ANALYSIS  ROUTINES  IN  WHICH  CAAE  THE  ANSWERS  DO',/, 

3 • NOT  REPRESENT  THE  TRUE  CHI-SQUARE  SOLUTUON') 

STOP 
END 


60 
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( JAN  73  » ^ ^ OS/360  FORTRAN  H 

COMPILER  OPTIONS  -\nAHE=  M A I N , 0PT*02 , LI NECNT  = 5 7 , S I ZE=0000K  , 

'SOURCE»EBCOICiNOLIST,OECK,LOAOf HAPiNOEOIT, lOtNOXREF 
SJBROUTINE  LSCPHI 
C 

C THIS  ROUTINE  CALCULATES  THE  THEORETICAL  FUNCTION,  RESIDUALS, 

C AND  CHI-SQUARE.  THE  USER  MUST  SUPPLYTHE  FUNCTION  PHIFNC 
C 

COMMON  A I 50,501 , C ( 50 1 , DPR  I 50 » , FMT I 5 ) , PR  I 501 , PS AVR ( 50 » , Z ( 50 » , P ( 50 1 
COMMON  PHI (10241 , RES  to (102^) ,X( 102A) ,Y(  10241 ,M(  10241 
COMMON  IFCOR , I F PLOT, I NOE Z , NB , NBADPT , NCHAN ,NE , NP , NPSTE P,  NRP , 

I NREAO,NWRI TE,NPUNCH 

COMMON  AA, CHI , CHIB, CHIA,0EGF 
C 

C CALCULATE  CHI, PHI, AND  RESIO 
C 

CHI=0.0 

00  10  I«l, NCHAN 
PHI  m=PHIFNC(  II 
R = SIO(  n = ( Y(  I l-PHI  ( I I l♦W(  I I 
CHI=CHIfRESIO( I I^^Z 
10  CONTINUE 
C 

C DIVIDE  CHI  BY  THE  # OF  DEGREES  OF  FREEDOM 
C 

CH  I=CHI/0EGF 

RETURN 

END 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  MA I N , 0PT=02 , L I NECNT=5 7 , S I ZE=OOOOK , 

SOURCE , EBCDIC »NOL 1ST, DECK, LOAD, MAP, NOEDIT, lOtNOKREF 
SUBROUTINE  ELTSI I I 


THIS  ROUTINE  CALCULATES  THE  DERIVATIVES  OF  PHIII). 
THE  USER  MUST  SUPPLY  THE  FUNCTION  DPHIII,J) 


COMMON  A(50,50) ,C( 50» ,DPRI50» ,FMT( 5) ,PR(50) , PS AVR 1 50 1 , Z ( 50  I , P ( 50 » 
COMMON  PHI  ( 132<»)  ,RESID(  132<*I  ,X(  13241  ,YI  1 324 1 , W ( 1 324 » 

COMMON  IFCOR, IFPLOT, I NDE Z , NB , NB AOPT , NCHAN , NE , NP , NPSTE P , N RP , 

I NREAD,NWRITE,NPUNCH 

COMMON  AA, CHI , CHIB, CH1A,0EGF 
00  10  J=NB,NE 
Z( JI=DPHI ( I, JI 
10  CONTINUE 
RETURN 
END 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  MA I N , 0PT=02 t L INECNT*57 , S I ZE=0000K, 

SOURCE  » EBCDIC  r NOLI  ST, DECK »LOAO» MAP »NOEO IT » 10,N0XREF 
SUBROUTINE  READIN 

♦♦♦  T^IS  ROUTINE  READS  IN  ALL  THE  INPUT  PARAMETERS,  DATA,  AND  SETS  UP 
INTERNAL  CONTROL  PARAMETERS. 

COMMON  A (50,501 ,C ( 50 » ,DPR I 501 ,FMT(5» ,PR(50) , PSAVR ( 50 1 , Z ( 50 1 , P ( 501 
COMMON  PHI ( 10241 ,RES ID ( 1024  I , X ( 1024 » , Y ( I 0241 , W( 1024  I 
COMMON  IFCOR, IFPLOT, INOEZ,NB,NBADPT,NCHAN,NE,NP,NPSTEP,NRP, 

I NREA0,NWRI TE,NPUNCH 

COMMON  AA,CHI ,CHI B,CHIA,DEGF 
DIMENSION  NBAD(  1 024 1 , YRUN ( 1024 ) ,TITLE(20  I 
EOU I VALENCE  (X( II , YRUNI 1 1 , NBAD ( L I I , I W ( L I , TI TLE ( 1 1) 
HRITE(NWRITE,LOOl) 

READ(NREAD,IITITLE 
WRirE(NWRITE,2l  TITLE 

***  READ  AND  WRITE  THE  CONTROL  PARAMETERS. 

RE A0(NREAD,3I  NCH AN.NBADPT , NRUNS , NRP , NPSTEP, I FPLOT , I FCOR 
WR ITEINWRITE, 10021 

WRITE(NWRITE,4I  NCHAN , NBADPT, NRUNS , NRP, NPSTEP , IFPLOT , IFCOR 

CALCULATE  THE  INTERNAL  CONTROL  PARAMETERS. 

IFINPSTEP.EO.CI  NPSTEP=l 
NPARAM=NRP 
NP*NRP 
NB  = l 
N£=NRP 
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C 

READ  ANO  WRITE  THE  PARAMETERS  FOR  THE  FUNCTION. 
C 

REAO(NREAO»5J (PRI II f I=l,NRPI 
W«ITEINWRITE»1003I 
WR  I TEINWRI  TE,6H  I t PR  (II » I^ltNRPI 
DO  40  1*1, NRP 
PI  I l*PR(I I 
40  CONTINUE 

C 

C + INITIALUE  FOR  SUMMING  NRUNS  OF  DATA. 

C 

DO  130  I=l,NCHAN 
Y(I  1 = 0.  0 
130  CONTINUE 

C***  GET  DATA  ANO  SUM  IT. 

DO  150  K=l, NRUNS 

READ(NREAD,FMT I (YRUN( I I ,I  = l,NCHANI 
DO  140  I=l,NCHAN 

Y( I I = YI I If YRUNI I I 
140  CONTINUE 

150  CONTINUE 

C 

C***  GET  THE  X AND  W VALUES  FOR  EACH  DATA  POINT 

PEAD(NREAD,FMTI  ( X ( I I , I * I , NST AR I 
Rr AD(NREAD,FMTI  ( WII I , I* I t NST AR I 

C 

C***  TEST  FOR  BAD  DATA  POINTS. 

C 

180  IF(NBADPTI210,210,190 

190  WRITE(NHRITE, 10061 
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c*** 


200 

c*<=* 

210 


READ  AND  WRITE  CHANNEL  « OF  BAD  POINTS. 
REA0(NREA0,9MNBAD(  1 ) , I-L»NBA0PT) 
W^ITE(NWRITE»9) (NBAD ( II » I» 1 .NB AOPT I 
DISCARD  BAD  DATA  POINTS. 

03  200  I=l,NBADPT 
J==NBA0(  1 I 
VI  Jl  = 0.  3 
W(J)»0.0 
CONTINUE 
HIITE  DATA. 

WR I TE (NWRITE ,1007  I 

WR I TEINWRI TE, 121 ( Y(I I ,I«I ,NCHAN) 

WRITEINWRITE, 13381 

WR I TE( NWRITE, 121  ( X( I I ,Isl,NCHANI 

WRl TEINWRITE, 10091 

WRl TEINWRI  TE  ,121  (W( I)  ,I-l,NCHANI 

WRITEINWRITE, lOOl I 


♦♦♦  CALCULATE  THE  H OF  DEGREES  OF  FREEDOM  FOR  THE  FIT. 


DEGF=NCHAN-NBAOPT-NRP 


***  CALCULATE  THE  INITIAL  THEORETICAL  FUNCTION  AND  CHI-SQUARED. 


CALL  LSQPHI 

chib=:hi 

***  SET  CHI  AFTER  AN  ITERATION  TO  100. 0 IN  CASE  PROGRAM  CAN  MAKE  NO 
SUCCESSFUL  ITERATIONS.. 

CHIA=100.0 

lOOl  FORMAT  IlHl,/,lH-» 

I FORMAT  I20A^» 
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2 FDRMAT  (20X,20AA» 

3 FDRMAT(7I5) 

1002  FDRMAT (IH-*  NCHAN  NBAOPT  NRUNS  NRP  NPSTEP 

I • I FCOR ' ) 

A FDRMAT  ( IX, /I  8» 

5 FDRMAT(8F10,4) 

1003  FDRMAT  ( IH- ' P AR AME T E R NUMBE R ' , 5X , ' P AR AME T ER • ) 

6 FDRMAT  ( I 4 X , I 3 , A X , F 1 0 . 2 ) 

7 FORMAT (4 ( I 5,F 10.5 ) I 

9 FDRMAT(14I5) 

1006  FDRMAT  ('-CHANNEL  NUMBER  OF  THE  BAD  DATA  POINTS') 

1007  FDRMAT  ( IH I , / , I H- , 30X , ' EX PER  I MENT AL  D AT  A ' , / / ) 

1008  FDRMAT( IHl ,30X, 'X  VALUES  FO  THE  EXPERIMENTAL  DATA') 

1009  FORM  AT ( IHl  , 30X , ' W VALUES  FO  THE  EXPERIMENTAL  DATA') 
12  FDRMATdOI  IX,E12.6)) 

RE  TURN 
END 


IFPLDT 
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COMPILER  OPTIONS  - NAME*  MA IN »0PT*0 2 , LI NECNT»57 , S I ZE»OOOOK , 

SOURCE,EBCOIC.NOL1ST,OECK>LOAO>MAP,NOEOIT» lOrNOXREF 
SUBROUTINE  TESTS 
C 

C***  THIS  ROUTINE  DOES  THE  TESTING  OF  THE  CHANGES  IN  THE  PARAMETERS 
C 

COMMON  A (50, 50) ,C ( 50 ) , OPR ( 50) , FMT ( 5 ) , PR  I 50) , PSAVR I 50 ) , Z I 50 ) , P ( 501 
COMMON  PHI  U024)  , RE S 1 0 ( 1 024  ) ,X  ( 1024 ) , V 1 1 024 ) , W ( 1024 ) 

COMMON  IFCOR, IFPLOT, INDEZ,NB,NBAOPT,NCHAN,NE,NP,NPSTEP,NRP, 

I NREAD.NWRI TEfNPUNCH 

COMMON  AA, CHI , CHIB, CHlAfOEGF 
DIMENSION  WE0(4) ,ZED(4) 

C***  AND  PERFORMS  THE  NECESSARY  CHANGES. 

C***  INITIALIZE  FOR  SUMMING, 

DT=0.0 

C***  SUM  AND  SAVE  THE  PARAMETERS. 

DO  10  I=NB,NE 

DT=0T+0PR ( I )*C( I ) 

PSAVR( I l=PR( I ) 

10  CONTINUE 

C***  TEST  FOR  DT  < ZERO. 

IF (DT) 20,40,40 

C*^*  CHANGE  THE  SIGN  OF  THE  DELTA  P»S  IF  DT<  ZERO. 

20  DO  30  I=NB,NE 

DPR(  I)=-OPR(I) 

30  CONTINUE 

DT=-DT 

C^**  USE  THE. FULL  LENGTH  OF  THE  CHANGES  VECTOR  TO  FIND  A NEW  CHI. 

40  AA=l.O 

C***  FIND  THE  NEW  CHI. 

CALL  CALS2 

C***  TEST  FOR  NEW  CHI  SMALLER  THAN  OLD,  IF  NOT  GO  TO  STATEMENT  90, 
IFICHIB.LE.CHI ) GO  TO  90 

C***  TEST  FOR  THE  CHANGE  IN  CHI  LESS  THAN  TERMINATION  VALUE. 
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50  IFnCHIB-CHD.GT.O.  IE-08J  GO  TO  80 

C***  REACHES  HERE  IF  TERMINATED, 

IMDEZ=0 

SET  CHI  AFTER  ITERATION  EQUAL  TO  CHI, 

60  CHIA=CHI 

C***  WRITE  OUT  CHANGES  IN  PARAMETERS#  LENGTH  OF  VECTOR, CHI  BEFORE, 

C***  AND  CHI  AFTER. 

WRI TEINWRITE ,l ) I I ,OPR( 1 1 , 1»NB,NE) 

WRITE(NWRITE,2) AA,CHIB,CHIA 
C***  SET  CHI  BEFORE  TO  CHI  AFTER. 

70  CHIB=CHIA 

RETURN 

C***  REACHES  HERE  IF  NOT  CONVERGED  YET  AND  GOES  BACK  TO  WRITE  OUT. 

80  INDEZ=-l 

GO  TO  60 

C***  REACHES  HERE  IF  FULL  LENGTH  OF  VECTOR  WILL  NOT  LOWER  CHI  AND 
C***  TRIES  THE  PARABOLAS  DESCRIBED  IN  WRITE  UP. 

90  S=CHI 

AA=0.5 
CALL  CALS2 
ZED(n=CHI 
WE'Om=AA 
N=l 
K*0 

AA=0T/(S-CHIB»2.0*OTI 

100  K*Kf1  

C***  IF  THE  FACTOR  MULTIPLYING  THE  CHANGES  VECTOR  IS  TOO  SMALL  SKIP  IT 
C***  FOR  THE  CHI  AFTER  CALCULATION. 

IF(AA,LT,0.1E-01I  GO  TO  llO 

CALL  CALS2 

N=N*1 
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ZED(N)>CH1 

WED(N)=AA 

ILO  IP(K-2I 120»130.153 
120  AA=l,0/0T 

60  TO  100 

130  ZE0l*CHIBi-S-2.0*ZED(ll 
lF(ZEDllL40tl50f L40 

UD  AA=(  Sf3.0*CHIB-4.0*ZEDm  l/(4.a*ZEDll 
GO  TO  100 
150  AA=0.5 

:hi«zeoc  1) 

00  160  I=«1»N 

IFlCHI.GT.ZEOUn  GO  TO  160 
CHI=ZE0(I) 

AA=HEO( I ) 

160  CONTINUE 

€♦♦♦  TEST  THE  SMALLEST  VALUE  GIVEN  BY  THE  PARABOLAS  TO  SEE  IF  IT  IS 
C***  SMALLER  THAN  CHI  BEFORE,  IF  NOT  RESET  PARAMETERS  AND  GO  TO  NEXT 
C***  SUBSPACE. 

IFICHIB.LE.CHU  GO  TO  170 
CALL  CALS2 
60  TO  50 

C***  . RESET  PARAMETERS. 

170  DO  100  I»NB,NE 

PR(I»=PSAVR( I » 

180  CONTINUE 

AA=0.0 
CALL  CALS2 
l'JDEZ=0 

190  WRITEINWRITEflOOl) 

WRITEINWRlTE.m  I ,OPR  1 1)  , I =NB , NE ) 
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1 

2 

I 00  I 


FJRMAT  0P( , I2,2HI  = ,E12.5,2<) ) 

FORMAT  AA=,Fl4.7f&X,bHCHll3=,El4.7,3X,5HCHIA=,El4.7,/l 

FDRMAT  ( *-rHF:  FOLLOUING  UPR(I»  PRODUCED  A DIVERGENT  STEP  WHICH 
I 'COULD  NOT  BE  FIXED  * I 

RETURN 
END 
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COMPILER  OPTIONS  - NAME*  MAI N» 0PT*02 t LI NECNT=57 i S I ZE'OOOOKt 

SOURCEi EBCDIC t NOLI ST t DECK tLOAOtMAPfNOEO IT » 10 tNOXREF 
SUBROUTINE  GRINOlINVl 
C 

C***  THIS  ROUTINE  SETS  UP  THE  MATRIX  OF  OERIVATIVESt  THE  CONSTANTS  IN 
C***  THE  NORMAL  EOUATIONSt  THEN  CALLS  THE  MATRIX  INVERTER,  AND  FINALLY 
C***  SHIFTS  THE  CONSTANTS  AND  THE  SOLUTION  TO  MATCH  THE  PARAMETERS. 

C 

COMMON  A(50,50) ,C( SO) ,OPR 1 50  I , FMT ( 5) , PRI 50) , PSAVR 1 50 ) , 21 50) ,PI50) 
CDMMONI  PHI  1102^) , RES  ID ( 1024)  ,X  11024)  ,Y  11024 ),  MI  1024) 

COMMON  IFCOR, IFPLOT, INDEZ,NB,NBAOPT,NCHAN,NE,NP,NPSTEP,NRP, 

I NREAO,NWRI TE,NPUNCH 

COMMON  AA tCHI  iCHIBfCHI A,OEGF 
DIMENSION  BI50) 

EQUIVALENCE  ( OPR  I 1 ) , B ( 1 ) ) 

C***  INITIALIZE  FOR  SUMMING. 

DO  20  1=1, NP 


DO  10  J=I,NP 

A1 I , J) *0.0 


10 

CONTINUE 

B( I )*0.0 

20 

CONTINUE 

SET 

UP  MATRIX  AND  CONSTANTS. 

DO' 

70  l=l,NCHAN 

c*** 

TEST  FOR  BAD  DATA  POINTS 

IF! Yl I) )40, 70,40 

40 

CALL  ELTSI I ) 

DO  60  J=1,NP  ' 

DC  50  K=J,NP 

A(  J,K)*A(J,K)f  ZUl^-ZIK) 
50  CONTINUE 

BIJ)=BI J)tZ(J)«RESIOII) 

60  CONTINUE 

70  CONTINUE 
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C***  SET  LDWFR  HALF  OF  MATRIX  (THE  MATRIX  IS  SYMETERIC). 
D3  93  J=I,NP 

[)0  8 0 K=J,NP 

A(K,JI=A( J,K> 

80  CONTINUE 

C SAVE  THE  CONSTANTS. 

C( J)=OI J) 

90  CONTINUE 

C**«  CALL  THE  MATRIX  INVERTER. 

CALL  INVMAT  (NP,INV,ZZZ) 

120  .RETURN 
END 
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COMPILER  OPTIONS  - NAME=  MA  I N , 0PT=02 f L I NECNT*5 7 » S I ZE=OOOOK, 

SOURCE t EBCDIC t NOLI  ST » DECK, LOAD, HAP »NOEO IT  I IO»NOXREF 
SUBROUTINE  INVMAT  ( N, I NV , P I VOT » 

COMMON  A I 50, 501 , C I 50 » , OPR C 50 1 ,FM T ( 5 I , PR  I 501 , PS AVR ( 50 1 , Z I 50 » , P ( 50 1 
COMMON  PHI  I I02<r»  ,RESIOIL02<r ) ,X  ( 1024)  ,Y(  1024)  ,W(  1024) 

COMMON  IFCOR,IFPLOT,INOEZ,NB,NBAOPT,NCHAN,NE,NP,NPSTEP,NRP, 

1 NREAOfNWRITEtNPUNCH 

COMMON  AA , CHI  , CHIB, CHIAtOEGF 

DIMENSION  B(50) ,IN0C(50) , I NOR ( 50 ) , I P I VOT (501 

EOU I VALENCE  ( JROM,IROW),  UCOLUM,  I COLUM I , ( AMAX,  SWAP,T) 

EQUIVALENCE  (OPR( I ) , B ( 1) ) 

C«‘4‘*  MATRIX  INVERSION  WITH  ACCOMPANYING  SOLUTION  OF  LINEAR  EQUATIONS. 


^♦♦♦♦♦♦♦♦note.  pivot  is 

C********NOTE.  INV=l  IF 


RETURNED  AS 
THE  INVERSE 


ZERO  IF  MATRIX 
IS  DESIRED. 


IS  SINGULAR. 


Z*** 


10 


20 


30 

40 

C*** 


INITIALIZATION 
DO  10  1=1, N 

IPI VGT(I 1=0 
CONTINUE 
DO  100  1 = 1, N 

SEARCH  FOR  PIVOT  ELEMENT. 

AMAX=0.0 
DO  40  J=1,N 

IF! IPIVOTI J) .EQ. I ) GO  TO  40 
DO  30  K=1,N 

IF( IPIVOTIK)-!  )20,  30,  130 

IF ( ABSI AMAX) .GT.ABSIAI J,K) ) ) 

IROW*J 

ICOLUM=K 

AMAX=A( J,K) 

CONTINUE 

CONTINUE 

IPIVOTI I COLUM) = IPIVOTI I COLUM) f I 
INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON 
IFI IROW.EQ.ICOLUM)  GO  TO  60 


GO  TO  30 


DIAGONAL. 
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50 

63 

C*<=* 

70 


80 

90 


DO  53  J=l,N 

SWAP^AdROW,  Jl 
A(  IROW*  JI^AdCOLUM,  J) 

A( ICOLUM, J>=SWAP 
CONTINUE 
$WAP=B( IROW) 

B ( IROWI^BI ICOLUH) 

B( ICOLUMI=SWAP 
P I V0T=A( ICOLUM, ICOLUMI 
TEST  FOR  SINGULAR  MATRIX. 

IFIPI VOT.EQ.O.OI  GO  TO  130 
DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT. 
DPIVOT=l.O/PI VOT 
Al ICOLUM, ICOLUMI=l.O 
DO  73  J=lfN 

A( ICOLUM, J)sA( ICOLUM* J I OOP I VOT 
CONTINUE 

BdCOLUM»  = BdCOLUM)*DPI  VOT 
REDUCE  NON-PIVOT  ROWS. 

DO  90  J=l,N 

IF( J. EG. ICOLUMI  GO  TO  90 
T=A( J,IC0LUMI 
AU,  ICOLUMI^o.O 
DO  83  K=l,N 

AI J,KI=A(J,KI-A( IC0LUM,KI*T 
CONTINUE 

B( JI=B(JI-B(ICOLUMI*T 
CONTINUE 

SET  INDEX  IF  INVERSE  IS  DESIRED. 
IFdNV.NE.  II  GO  TO  103 
INDRI I l-IROH 
INDCd  I^ICOLUM 
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100  CONTINUE 

INTERCHANGE  COLUMNS  IF  INVERSE  IS  DESIRED. 
IFIINV.NE.il  GO  TO  130 
DO  120  1*1, N 

J*N*l-I 

IFIINDRIJ)  .EQ.INDCI  JM  GO  TO  120 
JRQW*INDR( J) 

JC0LUM*IN0C(JI 
DO  110  K*1,N 


SWAP*A(K, JROW) 

AIK,  JROW)=MK,  JCOLUH) 
AIK, JCOLUM»=SWAP 

110 

CONTINUE 

120 

CONTINUE 

130 

RETURN 

END 
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COMPILER  OPTIONS  - MAME=  M A I N , OP T = 02 , L I NECMT*5 7 , S I ZE=OOOOK , 

SOURCE .EBCDIC  f MOL  I ST, OECKtLOAD, MAP, NOEO IT,  IO,NOXREF 
• SUBROUTINE  CALS2 

THIS  ROUTINE  IS  CALLED  BY  TESTS  TO  SET  UP  THE  PARAMETERS 
AFTER  THE  P6RAMETERS  HAVE  BEEN  CHANGED,  AND  CALCULATE  CHI  BY 
BY  CALL ING  LSCPHI 


C 

C** 

C** 

C 


10 

c 

c 

c 

60 


COMMON 

COMMON 

COMMON 


I 


AI50,50I ,CI50» ,0PR(50I ,FMT(5»,PRt50»,PSAVRI50),Z(50»,PI50I 
PHI  1 102‘»»  , RESIDI 1024  » ,X  ( 102AI  , YI  1024)  ,W(  1324) 

IFCOR, IFPLQT, I NDE Z ,NB , NBAOPT , NCHAN ,NE , NP,  NPSTE P, NRP , 


RESTRICTED  PARAMETERS. 


NREAO.NWR 

ITE 

common 

AA ,CHI ,CH 

IB, 

PINO  THE  CHANGED 

RE 

on  10 

I=NB,NE 

PR  I I ) = PS AV 

RI  I 

PI  I )=PR( I ) 

CONTINUE 

F IND  THE  CHI . 

CALL  LSQPHI 
RE  TURN 
f \m 
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COMPILER  OPTIONS  - NAME=  HA  I N,0PT*02 1 LI NECNT»5 7» S I ZE»OOOOK , 

SOURCE, EBCDIC f NOLI  ST, DECK,LOAO,MAP, NOEDIT, ID, NOXREF 
SUBROUTINE  SCRIBE 
C 

C***  THIS  ROUTINE  DOES  THE  WRITING  OUT, 

C 

COMMON  A I 50, 50) ,C ( 50 ) ,OPR( 50 ) , FMT I 5 ) , PR  I 50 ) , PSAVR ( 50 ) , Z ( 50 ) , P I 50 ) 
COMMON  PHK  1024)  , RES  ID  1 1024)  ,X(1024),Y(  1024 ),  W 1 1024) 

COMMON  IFCOR, IFPLOT, INOEZ iNB , NBAOPT , NCHAN,NE , NP, NPSTEP,NRP , 

I NREA0,NWRI TE,NPUNCH 

COMMON  AA,CHI ,CHIB,CHI A,DEGF 
DIMENSION  PH  INI  512) ,YN( 512) , XVEL ( 512 ) 

WRI TEINWRI TE, 1 ) 

WRITE(NWRITE,2) 1 1 ,PRI I) , I>1,NRP) 

20  WRITEINHRI TE,3) 

HRITE(NWRITE,4)  IPHKl)  ,Isl,NCHAN) 

1 FORMAT  I ‘-PARAMETER  NUMBER* , lOX, 'PARAMETER* ) 

2 FORMAT  ( 14X, I 3,5X,E14. 7) 

3 FORMAT  I 1H1,/,1H-,30X, 'THEORETICAL  FIT',//) 

4 F0RMATI8 12X,E13.7) ) 

30  RETURN 

END 
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COMPILER  OPTIONS  - NAME=  MA  I N ,OPT  = 02 » L I NECNT  = 57, S I ZE*OOOOK, 

SOURCE, EBCDIC, NOLI  ST, DECK, LOAD, MAP»NOEO IT f ID,NOXREF 
SUBROUTINE  ERRHAT 
C 

C***  THIS  ROUTINE  DOES  THE  ERROR  ANALYSIS. 

C 

COMMON  AI50,50; ,CI50I ,DPR(50» ,FMTI5I ,PR(50» ,PSAVR(50» ,ZI 50I,P(  501 
COMMON  PHI ( 1024) , RES  ID (1024) ,X I L024) , Y ( 1024) , W I 1024) 

COMMON  IFCOR,IFPLOT , INDEZ,NB,NBADPT, NCHAN,NE,NP,NPSTEP,NRP, 

I NREA0,NMRI TE,NPUNCH 

COMMON  AA, CHI , CHIB, CHIA,DEGF 

DIMENSION  AC  SR S I 50 ), AMS (50) ,F I NT ( 50 ) , S 1 6 I ( 50 ) 

C***  TEST  FOR  NO  SUCCESSFUL  ITERATIONS. 

IF (CHI  A. EQ. 103. 0)  GO  TO  163 
PIE=3. 1415927 

C**»  INITIALIZE  FOR  FINDING  THE  DATA  POINT  WITH  MAXIMUM  WEIGHTED 
C^**  RESIDUAL. 

TAX*0.0 

C***  FIND  CHANNEL  U OF  MAXIMUM  WEIGHTED  RESIDUAL  AND  WRITE  IT  OUT. 

00  10  I=l,NCHAN 

Q»ABS(RESI0( I ) ) 

IF(O.LT.TAX)  GO  TO  10 

ITRACK=i 

TAX=0 

10  CONTINUE 

WRITEINWRITE, DITRACK 

C***  FIND  THE  PROBABILITY  FOR  THE  FITTED  CHI. 

AR=SORT( 2.0*CHIA*DEGF)-S0RTI2.3*DEGF-1.3) 

PROB=ERFC( AR)/2.0 
NCH=NCHAN-NBADPT 

C***  WRITE  OUT  STATISTICAL  INFORMATION. 

WR ITEINWRITE ,2)NCH,NRP,CHI A,PROB 
C*-**  TEST  CHI  FOR  LESS  THAN  l.O  AND  SET  TO  ONE  IS  SATISFIED. 

IF (CHI A.LT. l.O)  CHIA=1.3 
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C***  CALL  GRIND  TO  FIND  THE  ERROR  MATRIX. 

20  CALL  GRINOm 

C***  FIND  THE  STANDARD  DEVIATIONS  ON  THE  FITTED  PARAMETERS. 
D3  30  I=l,NRP 

AC  SR  SI I)=SQRT(A(I,i I*CHIA» 

30  CONTINUE 

C***  FIND  THE  CORRELATION  COEFFICIENTS. 

II=NRP-l 
DO  50  1=1, II 
JJ=I^l 

SIG=ACSRS(  I I 
DO  40  J=JJ,NRP 

A(I,JI=At  I,  JI*CHIA/(SIG*ACSRSU»  I 
40  CONTINUE 

50  CONTINUE 

C***  WRITE  OUT  THE  PARAMETERS  AND  THEIR  STANDARD  DEVIATIONS. 
WRI TEINWRITE ,31 

WRI TEI NWRI TE,4) ( I ,PR( I 1 ,ACSRS( II  , l = l,NRPI 
C***  TEST  FOR  CORRELATION  COEFFICIENTS  DESIRED. 

100  IF  ( IFCOR)  no,  150, 120 

C***  IF  IFCOR  < ZERO  PUNCH  OUT  PARAMETERS, 
no  WRITE(NPUNCH,Ill ( PR ( I 1 , ACSRS II 1 , I = l,NRPI 
120  WRITE(NWRITE,12I 

C***  WRITE  OUT  CORRELATION  COEFFICIENTS. 

I I =NRP-l 
DO  140  1 = 1, II 
jj=m 

WRITE(NWRITE,L31(I,J,A(I,J1, J=JJ,NRP) 

C***  TEST  FOR  PUNCHING. 

IF( IFC0R)130,140,140 

130  WRIT£(NPUNCH,14)(A(I,JI , J=JJ,NRP) 
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CONTINUE 


a 


140 

150 

160 

1 

2 

3 


11 
12 
13 
1 4 
15 


RETURN 


k^RITEINWKI  TE  ,151 

FORMAT  ('-CHANNEL  NUMBER  OF  THE  WORST  DATA  POINT  IS  ',151 
FORMAT  I '-CHI-SQUARED  FOR', 15,'  CHANNELS  AND', 15,'  PARAMETERS  IS', 
1 Fll.7,/,'  THE  PROBABILITY  FOR  WHICH  IS',E16.7) 

FORMAT  ( IHl, /, '-PARAMETER  NUMBER ', lOX ,' PARAMETER ', lOX ,' STANDARD  ', 
1 'DEVIATION'! 

FORMAT  (14X,I3,5X,E14.7,14X,E14.7) 

FORMAT  ( 1X,E14.7,15X,E14.7! 

FORMAT  (5E16.7) 

FORMAT  ( IHl, /, IH-,30X, 'CORRELATION  COEF F I CENT S ' I 
FORMAT  (1H-,6(3H  A ( , I 2 , 1 H , , I 2 , 2H 1 = , F 1 0 . 7 I ) 

FORMAT  (8F10.51 

FORMAT  ('  PROGRAM  COULD  MAKE  NO  SUCCESSFUL  ITERATIONS'! 

RE  TURN 


END 
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0S/36D  FORTRAN  H 


COMPILER  OPTIONS  - NAME*  MA IN,0PT»02, LIN£CNT=57tS I ZE=OOOOKt 

SOURCE, EBCDICt NOLIST, DECK, LOAD, MAP, NOEDIT, ID,NOXREF 
SUBROUTINE  PLOTS 


♦♦♦  THIS  ROUTINE  DOES  THE  PLOTTING  OF  THE  DATA  AND  THEORETICAL  FIT, 


VC' 


COMMON  A (50, 501 ,C(501 , DPR ( 50 ) , FMT ( 5 I , PR ( 50 » , PS AVR 1 50 » , Z ( 50  I , P ( 501 
COMMON  PHI ( 10241 , RES  10(1024  I ,X(  1024) ,Y( 10241 ,W(  1024) 

COMMON  IFCOR, IFPLOT, INOEZ,NB,NBAOPT,NCHAN,NE,NP,NPSTEP,NRP, 

I NREAO,NHRI TE,NPUNCH 

COMMON  AA, CHI , CHIB, CHIA,OEGF 
DIMENSION  ARR(126) 

WR ITEINWRI TE,1 ) 

DATA  BLANK, STAR, ZERO, PLUS, BAO/IH  , IH*,  IHO,  IH4-,  LHB/ 

C**+  INITIALIZE  FOR  FINDING  MAXIMUM  AND  MINIMUM, 

AMAX=PHI I 1) 

AMIN=AMAX  . 

C***  INITIALIZE  THE  PLOTTING  ARRAY  TO  BLANKS, 

DO  10  1*1,126 

A^R( I )*BLANK 
10  CONTINUE 

£♦♦♦  FIND  MAXIMUM  AND  MINIMUM  OF  EITHER  THE  DATA  OR  THE  FIT. 

DO  30  I=l,NCHAN,NPSTEP 


IF (Y( I ),EQ.O.O)  GO 

TO  20 

IF( Y( I ) .GT.AMAX) 

AMAX=Y( 1 ) 

IFI Y( I 1 .LT.AMIN) 

AMIN*Y( I ) 

20 

IF(PHKI). GT.AMAX) 

AMAX=PHI (I ) 

IF(PHI( I ) .LT.AMIN) 

AMIN* PHI ( I ) 

30 

CONTINUE 

c*** 

FIND  THE  BIN  SIZE. 

DB IN=( AMAX-AMIN) /124.0 

c*** 

FIND  WHERE  TO  PUT  THE  CURVES. 

DO  80  l=l,NCHAN,NPSTEP 

II=IFIX( (PHI (I)-AMIN)/0BIN)»2 
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IFtYd  ).NE.O.O»  GO  TO  43 
J J = l 

&RR( I »=BAO 
GO  TO  6C 

JJ  = IFIX((Y(l  >-AMlNI/OBlN)*-2 

IF ( I I .NE. J J)  GO  TO  53 

ARRII I >=PLUS 

GO  TO  70 

ARRI JJ)=ZERO 

ARR(m=STAR 

WRI TEINKRI TE  ,2>  I , ( ARR(K) ,K=l, 126) 

ARR  ( I I ) = BLANK 
ARR  ( J Jl  =BLANK 
CONTINUE 
WR I TEINWRI TE  ,3) 

FORMAT  ( IHl) 

FORMAT  ( IXt I 4, 2X,  126A  I ) 

FORMAT  ( *-0=EXPERIMENTAL  POINT,  ♦ = FITTEO  POINT,  ♦•=  COI NC I DENT  ', 

1 'EXPERIMENTAL  AND  FITTED  POINTS,  B=BAO  EXPERIMENTAL  ', 

2 'POINT') 

RE  TURN 

END 
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OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME*  MA IN, OP T=02 » LI NECNT*57 t S I Z E=OOOOK t 

SOURCE, EBCDIC ,NOL I ST ,OECK,LOAO, MAP, NOED I T , IO,NOXREF 
FUNCTION  PHIFNCm 

COMMON  AtSO.bCI ,C150) , OPR (501 ,FM T 1 51 ,PR ( 50) , PSAVR 150) , Z ( 50 ) ,PI 501 
COMMON  PHI (1024) , RES  ID (1024) ,X 1 1024) ,Y( 1024) ,W(  1024) 

COMMON  IFCCR , I F PLOT , I NOE Z , NB , NBAOPT , NCHAN , NE , NP , NPSTEP , NRP , 

1 NREAO,NWRI TE,NPUNCH 

COMMON  AA, CHI , CHIB, CH1A,DEGF 
DIMENSION  PSXI50) 

C 

C THIS  SUBROUTINE  CALCULATES  THE  THEORETICAL  FUNCTION  FOR  THE  FITTING 
C PROCESS  . SPECIFICALLY  PHI  «SUN  ( P (K ) *EX  P( -P  (K«- 1 )*X  ( I ) )*COS  ( P ( K »2  )4X(  I ) ) 
C WHERE  THE  SUM  IS  OYER  THE  TOTAL  NUMBER  OF  FUNCTIONS 
C 

XI*X( I ) 

DO  10  K=1,NP 
PSX(K)=P(K)*XI 
10  CONTINUE 

PH  IFNC=0.0 
00  20  J=l,NP,3 

20  PHIFNC*PHIFNC^P( J)*EXP(-PSX(J*  1) ) *COS I PS X (J ♦ 2 ) ) 

RETURN 

FNO 


APPENDIX  B 


1 152 


t I 


( JAN  73  I 


OS/360  FORTRAN  H 


COMPILER  OPTIONS  - NAME=  MA I N, OPT =0 2 , L I NECNT *5 7 , S I ZE*OOOOK f 

SOURCE tEBCOICf NOLI  ST » DECK, LOAO»MAP,NOEOIT» ID>NOXREP 
FUNCTION  DPHK  I , J) 

COMMON  A (50, 50) ,C ( 50 ) ,DPR ( 50 1 , FMT ( 5 1 , PR ( 50 ) , PSAVR ( 50 1 , Z ( 50 ) ,P ( 50) 
COMMON  PHI  I 1024)  ,RESID(ID24)  ,XU024)  , Y(1024)  ,M(1024) 

COMMON  IFCOR, IFPLOT, INDEZ,NB,NBAOPT,NCHAN,NE,NP,NPSTEP,NRP, 

I NREADfNWR I TE,NPUNCH 

COMMON  AAtCHI fCHlBfCHI AfOEGF 
DIMENSION  PSX(50) 

DIMENSION  SAVS(50) ,SAVC(50) 


C 

C THIS  ROUTINE  CALCULATES  THE  DERIVATIVES  OF  THE  THEORETICAL  FUNCTION 
C WHERE  THE  FUNCTION  IS  GIVEN  BY 

C PHI  * SUM(  PI  K)*EXP(-P(KU  )♦<(  I ) )*C0S(P(K«^2)*X(  I ) ) ) WHERE  THE 
C SUM  IS  OVER  ALL  FUNCTIONS 
C 


IF(J.NE.I)  go  TO  20 
XI=X(  I ) 

DO  10  K=l,NP 
PSX(K)=P(K)*XI 


10  CONTINUE 

K*0 

DO  15  L=l,NP,3 
K*Kfl 

CON=P(  L)*EXP(-PSX(L«^l)  ) 
SAVCIX )=C0N*CCS(PSXILf2) ) 
SAVS(K)*C0N»SIN(PSX(L*-2)  ) 
15  CONTINUE 

20  ITE=J-l 

K=MOO( ITE,3) f 1 
L = I J-l  )/3fl 
GO  TO  (l,2,3),K 
I OPHI=SAVC(L) /PI J) 

RF  TURN 
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DPHI=-XI ♦SAVC IL> 
RFTURS 

0PHI=-XI*SAVSIU 

R£  TURN 
END 
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